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ABSTRACT 


We conduct a rigorous examination of the nearby red supergiant Betelgeuse by drawing on the 
synthesis of new observational data and three different modeling techniques. Our observational results 
include the release of new, processed photometric measurements collected with the space-based SMEI 
instrument prior to Betelgeuse’s recent, unprecedented dimming event. We detect the first radial 
overtone in the photometric data and report a period of 185 + 13.5d. 

Our theoretical predictions include self-consistent results from multi-timescale evolutionary, oscil- 
latory, and hydrodynamic simulations conducted with the Modules for Experiments in Stellar Astro- 
physics (MESA) software suite. Significant outcomes of our modeling efforts include a precise prediction 
for the star’s radius: up Ro. In concert with additional constraints, this allows us to derive a new, 
independent distance estimate of 168*27 pe and a parallax of 7 = 5.95+0:58 mas, in good agreement 
with Hipparcos but less so with recent radio measurements. 

Seismic results from both perturbed hydrostatic and evolving hydrodynamic simulations constrain 
the period and driving mechanisms of Betelgeuse’s dominant periodicities in new ways. Our analy- 
ses converge to the conclusion that Betelgeuse’s ~ 400 day period is the result of pulsation in the 
fundamental mode, driven by the &-mechanism. Grid-based hydrodynamic modeling reveals that the 
behavior of the oscillating envelope is mass-dependent, and likewise suggests that the non-linear pul- 
sation excitation time could serve as a mass constraint. 

Our results place œ Ori definitively in the early core helium-burning phase of the red supergiant 
branch. We report a present-day mass of 16.5-19 M—slightly lower than typical literature values. 


Keywords: stellar evolution — red giants — stellar oscillations — numerical techniques 


1. INTRODUCTION 


Since November of 2019, the red 
a Orionis—popularly known as Betelgeuse—has experi- 
enced an unprecedented brightness drop of nearly 2 mag- 
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nitudes in the V band. The severity of this decrease and 
the deviation from its typical pattern of variability have 
sparked much public speculation about the physics re- 
sponsible and its likelihood of undergoing a cataclysmic 
event. 

To investigate these questions first requires an under- 
standing of the short-timescale behavior of variable red 
giants. Such stars are known to exhibit a complex spec- 
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trum of variability, where cyclic variations with differ- 
ent driving mechanisms occur over a range of timescales. 
Though we can explain and fully capture some pulsation 
physics in 1D stellar models (e.g., pressure and gravity 
modes; see review by Aerts 2019), other mechanisms 
are not well understood (Wood et al. 2004; Nicholls 
et al. 2009a). In this latter class fall many of the varia- 
tions we observe on human timescales, as such behavior 
is, with rare exception, too rapid to be explained by 
classical stellar evolution (Molnar et al. 2019). Mod- 
eling such processes may require 3 dimensions, time- 
dependent convection, or otherwise more sophisticated 
physical formalisms that are beyond the scope of typical 
1D stellar evolution programs. Nevertheless, 1D stellar 
models are among the most powerful devices for gain- 
ing insight on the sub-surface physics responsible for 
observed changes in real stars (Demarque et al. 2004; 
Pietrinferni et al. 2004; VandenBerg et al. 2006; Cordier 
et al. 2007; Weiss & Schlattl 2008; Dotter et al. 2008; 
Townsend & Teitler 2013; Paxton et al. 2018 and oth- 
ers). When conducted on a range of timescales, their 
calculations can be exploited to great effect. 

In red supergiants, the &-mechanism drives radial pul- 
sations in the hydrogen ionization zone, and simula- 
tions show the emergence of periods and growth rates 
of the dominant fundamental pulsation mode—typically 
on the order of years—both in linear and non-linear 
models, as shown in e.g. Li & Gong (1994), Heger et al. 
(1997), Yoon & Cantiello (2010), and Paxton et al. 
(2013). In addition to these, previous modeling work 
on a Ori and similar red supergiants (RSGs) includes 
Dolan et al. (2016), Wheeler et al. (2017), Nance et al. 
(2018), and Goldberg et al. (2020). 

In both Yoon & Cantiello (2010) and Paxton et al. 
(2013), models of rotating and non-rotating RSGs with 
approximately solar metallicity and initial masses of 
25M were found to exhibit pulsations on the order 1-8 
years. Obtaining frequencies of this magnitude required 
lowering the evolutionary timestep to a fraction of a year 
during helium burning. The limiting factor on these cal- 
culations was the emergence of supersonic radial veloci- 
ties in the envelope (see Section 6.6 in Paxton et al. 2013 
for more details on their example). 

A rigorous estimation of the model-derived fundamen- 
tal parameters of œ Ori was undertaken by Dolan et al. 
(2016). In particular, their models find a best estimate 
of 20*3 Mo for the progenitor mass. They also attempt 
to model the pulsation properties of o Ori, but find 
they were unable to reproduce the fundamental mode 
(FM) and first overtone (O1) frequencies with adiabatic 
models alone. They suggest that interplay between non- 
adiabatic pulsations and convection could be responsi- 


ble for some variability, noting that 3D simulations of 
similar red supergiants show the development of large- 
scale granular convection that can itself drive pulsation 
(Xiong et al. 1998; Jacobs et al. 1998; Freytag et al. 
2002; Chiavassa et al. 2011; Freytag & Chiavassa 2013; 
Dolan et al. 2016). 

Recent 1D modeling efforts in “The Betelgeuse 
Project" series and other works suggest that a past 
merger may be required to explain the present-day sur- 
face rotation of o Ori, which is more rapid than standard 
stellar evolutionary calculations including rotation can 
reproduce (Wheeler et al. 2017; Nance et al. 2018; Chat- 
zopoulos et al. 2020). The Nance et al. (2018) study also 
examines the star seismically, but the authors are pri- 
marily focused on rapid waves in the convection zone 
that might precede a cataclysmic event. This concept 
was also addressed in depth by Goldberg et al. (2020), 
who modeled the observable features of supernova events 
as a function of the point of onset during the stellar pul- 
sation. 

In this paper, we use a range of tools to investigate 
the variability of œ Orionis. We use the Modules for 
Experiments in Stellar Astrophysics (MESA, Paxton 
et al. 2011, 2013, 2015, 2018, 2019) stellar evolution soft- 
ware suite to generate both classical evolutionary tracks 
and short timescale, hydrodynamic simulations of stars. 
We likewise use the GYRE pulsation program to con- 
struct complementary predictions of the pressure mode 
(p-mode) oscillations in models of red giants (Townsend 
& Teitler 2013). 

'This paper proceeds as follows: In Section 2, we dis- 
cuss the current knowledge of a Ori's classical con- 
straints, including pulsation periods, evolutionary stage, 
radius, temperature, and distance. We present a 
lightcurve highlighting o Ori's recent behavior, con- 
structed from data collected from the American Asso- 
ciation of Variable Star Observers (AAVSO) and newly 
processed space-based photometry from the Solar Mass 
Ejection Imager instrument. In Sections 3, 4, and 5, 
we discuss our evolutionary, seismic, and hydrodynamic 
models, respectively. Section 6 concludes our analysis 
and presents best estimates of its fundamental param- 
eters based on detailed photometric analysis and com- 
prehensive, multi-timescale simulation. 


2. OBSERVATIONAL CONSTRAINTS 


a Ori is well studied interferometrically; together with 
R Dor and IRC 10216, it is among the stars with the 
largest angular diameters ever measured (Bedding et al. 
1997; Menten et al. 2012; Stewart et al. 2016). In their 
Table 3, Dolan et al. (2016) provide a clear summary of 
previous measurements. 
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Table 1. Processed SMEI photometry of œ Ori. Observations 
were corrected for systematics and averaged into 1-day bins. 
Errors calculated as simple shot noise. V mag is the same 
light curve, scaled to existing V-band data. The full data set 
is available in the online journal 


BJD- SMEI SMEI V V 
2400000 (d) mag error mag error 
52677.959995 0.3759 0.0037 0.5168 0.0037 
52678.983194 0.3849 0.0039 0.5330 0.0039 
52680.041678 0.3717 0.0043 0.5094 0.0043 
52680.959028 0.3801 0.0039 0.5244 0.0039 
52681.911649 0.3869 0.0035 0.5365 0.0035 
52682.934838 0.3908 0.0036 0.5436 0.0036 
52683.887465 0.3991 0.0035 0.5586 0.0035 
52684.875377 0.4032 0.0035 0.5662 0.0035 
52685.863269 0.4030 0.0035 0.5657 0.0035 
52686.851169 0.4068 0.0035 0.5727 0.0035 
52687.415625 0.4177 0.0093 0.5929 0.0094 
52689.038675 0.4142 0.0042 0.5863 0.0042 


The earliest interferometric measurement from 
Michelson & Pease (1921) resulted in an angular di- 
ameter of 47+5 mas at visible wavelengths, assuming a 
uniformly illuminated disk model. In recent years, it was 
realized that there were elevated layers of molecules and 
dust above the photosphere (e.g., Perrin et al. 2004), 
complicating the interpretation of diameter measure- 
ments. In the context of the recent dimming event of 
a Ori, Haubois et al. (2019) solved for the photospheric 
diameter, dust shell diameter and optical depth. They 
found a Uniform Disk diameter of 44.0+0.5 mas in their 
1.04 wm, bandpass, which had relatively little influence 
from molecular bands. This would be equivalent to a 
limb-darkened diameter of 46.0+0.6 mas using a linear 
limb darkening coefficient of ~0.5 (Claret & Bloemen 
2011; Hanbury Brown et al. 1974). However, the rel- 
atively simple dust model consisting of of a 64.7 mas 
diameter thin shell scattering 4.4 96 of the light means 
that the statistical error from that work is not fully 
representative of the model uncertainty. 

We adopt as observational reference the diameter 
from the recent work of Montargés et al. (2014) of 
42.28+0.43mas for the limb-darkened (i.e., physical 
photospheric) diameter. Those authors resolved the 
photosphere significantly past the first null in the visibil- 
ity curve, so were insensitive to low optical depth shells, 
unlike many of the other measurements. Additionally, 
the relatively high spectral resolution observations in the 
K band, away from main molecular absorption features, 


mean that this measurement is relatively unaffected by 
apparent molecular shells. 

Radius estimates are further complicated by uncer- 
tain parallax measurements, which are made difficult by 
variability and known >2au-scale asymmetries on the 
surface of the star both at optical and radio wavelengths 
(Young et al. 2000; Kervella et al. 2018). The revised 
Hipparcos astrometric solution gave an optical-only dis- 
tance of 15377? pc (van Leeuwen 2007). Combination 
of the Hipparcos data with radio observations captured 
by the Very Large Array (VLA) extended that distance 
out to 197+ 45 pc, which was also used by Dolan et al. 
(2016). 

The revised Hipparcos-only value is inconsistent at 
the 1.7c level with the most recent radio measurement 
of 200 pc (Harper et al. 2017), which took into ac- 
count both VLA and Atacama Large Millimeter Array 
(ALMA) observations but which was also significantly 
affected by “cosmic noise".! The star is well beyond the 
established brightness limit of Gaia, and data enabling 
a future parallax measurement were not collected in the 
first years of the mission. A parallax estimate of Betel- 
geuse is therefore not included in Gaia Data Release 2 
(Gaia Collaboration et al. 2016; Sahlmann et al. 2018). 
Given the very long time-baselines needed to overcome 
the effects of photospheric motions and variability, there 
is unlikely to be a reliable direct parallax measurement 
of Betelgeuse with « 1096 uncertainty in the near term. 

Estimates of Betelgeuse's mass are derived from mod- 
els and range from roughly 15-25Mo, with previous 
modeling work suggesting that a Ori is in the midst 
of its core helium-burning giant branch phase (Neilson 
et al. 2011; Dolan et al. 2016; Wheeler et al. 2017; Nance 
et al. 2018). However, while Dolan et al. (2016) state 
that its mass loss rate—the primary piece of evidence 
supporting the claim that it is on the red supergiant 
branch (RSB)—is “consistent with having recently be- 
gun core helium burning," they also note that a previous 
interaction of Betelgeuse with a binary companion could 
account for similar mass loss rates without necessitating 
that Betelgeuse currently exist on the RSB. Since nearly 
half of ~ 20Mo stars have a companion close enough to 
induce mass loss, this is, in fact, ambiguous (de Mink 
et al. 2014). It is demonstrated by Wheeler et al. (2017) 
that rotating models of a Ori do not produce reasonable 


1 “Cosmic noise” is an umbrella term used to describe the elevated 
dispersion of the residuals of the astrometric solution compared 
to the formal errors. It can include various physical effects such 
as source size, unresolved companions, unresolved properties of 
stars in the stellar models used for fitting, variability of the stellar 
parameters, and instrumental effects such as excess noise due to 
saturation. 
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Table 2. Observational best values, estimated ranges, and model-derived constraints (where indicated) for o, Ori. The tempera- 
ture constraints reflect the spectroscopically derived temperature from o Ori at its brightness minimum, which is not necessarily 
reflective of its mean temperature. However, even Levesque & Massey (2020)’s 100 K error bars accounting for decadal variations 
are more restrictive than the theoretical uncertainty imposed by modeled variations in amir. Though we quote a “best” radius 
and reference a wide range of values, in practice we do not impose any constraints on the radius when modeling. The range of 
possible radii derived from the models is smaller than the uncertainties reported by many observers. We quote the initial and 
present-day mass ranges preferred by our seismic models. Masses considered in the initial grid range from 10 to 26 Mo. 


Property Value Source Comment 

Tog 3600 + 25 K Levesque & Massey (2020) range extended by Gtheory to +200 K 
Angular Diameter 42.28+0.43 Montargés et al. (2014) Limb-darkened 

Radius upper limit ~ 1100Ro Dolan et al. (2016) data collated from many sources 
Radius lower limit 500Re Dolan et al. (2016) data collated from many sources 
Distance 197 + 45 pc Harper et al. (2008) parallax data adopted by Dolan et al. (2016) 
Period of variability 388 + 30 days Kiss et al. (2006) dominant, higher frequency; likely FM 
Period of variability 2050 + 460 days Kiss et al. (2006) lower frequency; likely LSP 

Period of variability (FM) 416 + 24 days this work SMEI+ V data; mode ID from GYRE 
Period of variability (O1) 185 + 14 days this work SMEI+ V data; mode ID from GYRE 
Period ratio O1/FM 0.445 + 0.041 this work SMEI+ V data 

Radius 764 + 116, —62Rọo this work 30 range; seismic analysis 

Initial Mass 18-21Mo this work median range; seismic analysis 
Present-day mass 16.5-19Mo this work median range; seismic analysis 
Distance 168.1 + 27.5, —14.9 pc. this work 30 range; seismic analysis 

Parallax 5.95 + 0.58, —0.85 mas this work 30 range; seismic analysis 


evolutionary predictions (a finding consistent with our 
present work), but they do not draw any specific con- 
clusion about whether the star is core helium burning. 
As it is impossible to measure either mass or evolu- 
tionary status directly, and the evidence regarding its 
phase is not definitive, we do not assume a particular 
evolutionary phase a priori in our models. Instead, we 
consider the relative probabilities that o Ori is in a par- 
ticular evolutionary stage based on (1) the masses of 
tracks that match the other observational constraints 
and (2) the duration of the possible evolutionary stages. 
The first-order, theoretical constraints on its mass 
and age are provided by the linear pulsation calcula- 
tions, which rule out any model in an evolutionary stage 
earlier than the RSB. From an observational perspec- 
tive, we note that Betelgeuse is far in the foreground 
of the known «10 Myr age young associations in Orion 
(GroBschedl et al. 2018), and it is not known to have 
kinematics consistent with ejection. In particular, its 
radial velocity of --21.94-0.5 km s^! is consistent with 
the ~+23 km s^! of typical high mass stars in the Orion 
OBI association (Morrell & Levato 1991; Famaey et al. 
2005), but would differ by —20 kms! if it had trav- 
elled 200 pc in ~20 Myr. The (U, V, W) space motion of 
Betelgeuse is (—22, —10,12) kms"! with respect to the 
Sun, which is (—11,2,19) kms~! with respect to the lo- 


cal standard of rest (Famaey et al. 2005; Schónrich et al. 
2010). The high W velocity in particular is of note, as 
it is discrepant at 3o from the kinematics of the young 
disk (Robin et al. 2003). If this high W velocity were 
due to ejection from a young association lying on the 
Galactic disk, now falling back through the disk due to 
vertical epicyclic motion, this would imply an origin of 
~50 Myr ago. With these proper motion estimates in 
mind, we are left with a few possible scenarios of vary- 
ing likelihood: (1) Betelgeuse was formed very recently 
in a region where there is no star formation; (2) it is 
250 Myr old, or (3) it underwent some kind of binary 
interaction that propelled its trajectory. Scenario (1) 
is not reasonable, and scenario (2) would be consistent 
with a mass below 10Mo—a possibility that is readily 
ruled out by our other constraints. We are thus left with 
the third scenario, which is likewise supported by obser- 
vations of Betelgeuse's present-day surface rotation and 
the inability of 1D, rotating models to reproduce it (see 
subsequent discussion). 

We construct an age-prior function that performs a 
Monte Carlo interpolation over a grid of stellar tracks 
with masses ranging from 16-26Mo (other parameters 
fixed; amur = 2.1) and a power law IMF. For two sets 
of realizations, we adopt a minimum age constraint of 8 
Myr and permit radii between 600 and 9005. In the 
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first statistical run, masses are heavily skewed towards 
the head of the distribution, peaking at 16Mo, and the 
bulk of the trials fall from 16-18M.. This indicates 
that the lower-mass regime is strongly statistically pre- 
ferred, which is consistent with our understanding of the 
prevalence of high-mass stars in general. In the second 
statistical run, the distribution peaks a bit higher, at 
18Mo, and tapers off rapidly beyond 17.5 and 19.5 in 
either direction. The number of trials that do not fall 
somewhere on the core helium-burning RSG is negligi- 
ble regardless of mass, though this is even more strongly 
the case for trials with masses between 17 and 19Mo. 

As we will conduct estimates of the stellar mass, 
and many other parameters, in several ways throughout 
this analysis, we treat the above statistical experiment 
strictly as sufficient evidence to assume that Betelgeuse 
is core helium burning in subsequent modeling. 

Recent spectral analysis of Betelgeuse presents an ef- 
fective temperature of 3600 + 25 K (e.g., Levesque & 
Massey 2020; Guinan & Wasatonic 2020). Levesque 
& Massey (2020) write that the difference between 
the spectroscopically-derived temperature measured in 
2004-5 and that measured during Betelgeuse’s bright- 
ness minimum in 2020 is at most a decrease of 100K, and 
at minimum, negligible. We note, however, that there is 
some disagreement on the reliability of the method by 
which Levesque & Massey (2020) derive their tempera- 
ture, with, e.g., Ireland et al. (2008) and Davies et al. 
(2013) suggesting it may be underestimated. We discuss 
this in more detail in Section 3. Adopting the results of 
Levesque & Massey (2020) essentially rules out convec- 
tive turnover as an explanation for its recent dimming, 
but surface temperature is less informative on other os- 
cillation driving mechanisms. 

Critically, the brightness of Betelgeuse varies in a sys- 
tematic way on at least two different timescales, and 
these periodicities were measured with good precision 
by Kiss et al. (2006) (and later corroborated by Chatys 
et al. 2019). The shorter occurs with a period of ~ 388 
days and the longer with a period of — 5.6 years (2050 
d). The period-luminosity relation depicted in Figure 6 
of Kiss et al. (2006) provides some evidence that the 388 
d brightness variation is caused by p-mode pulsation in 
the fundamental mode. This is likewise supported by 
various observational and theoretical considerations, in- 
cluding the position of the star on the log P-Mx dia- 
gram, where the absolute K brightness provides the ob- 
servational proxy for the stellar luminosity. T Kiss et al. 
(2006) also found that the shorter periods fit the theoret- 
ical calculations of Guo & Li (2002), forming an exten- 
sion to sequences B and C of the supergiant variables 
observed in the Magellanic Clouds that correspond to 


the FM and the O1 pulsation modes, respectively (Wood 
et al. 1999; Kiss & Bedding 2003; Soszynski et al. 2007). 
This also suggests that these variations correspond to 
p-mode pulsation. 

The longer, ~ 2050 d periodicity likely falls in a class 
of signal known as *Long Secondary Periods," or LSPs. 
These have been observed in multiple semiregular and 
red supergiant variables, but the cause of the LSP is 
still debated (Wood 2000; Chatys et al. 2019). Pro- 
posed mechanisms include rotational modulation caused 
by spots or a nearby companion followed by a dust 
cloud, among other possibilities (Wood 2000; Soszyński 
& Udalski 2014). Such signals were observed in the LMC 
supergiant population as “sequence D,” and the long pe- 
riods found by Kiss et al. (2006) extend that sequence to 
higher luminosities (Derekas et al. 2006). Among other 
things, rotational modulation was proposed as a possible 
mechanism for the LSP (Percy & Deibert 2016). How- 
ever, the rotational period of œ Ori has recently been 
estimated at Pot = 3148 yr, which is considerably 
longer than the LSP of the star (Kervella et al. 2018). 
Models in this work shed more light on the questions of 
mode classification and driving mechanism. 


2.1. Photometric Observations 


Both the ~ 400 day and 5.6 yr (2050 d) periods are 
visible in Figure 1, which shows the longitudinal bright- 
ness variations of œ Ori over the last 90 years. These 
visual brightness estimates were collected in large part 
by amateur observers and archived by the American As- 
sociation of Variable Star Observers (AAVSO). 

Examining Figure 1 more closely, we see that the am- 
plitude of the brightness drops corresponding to the 
~ 400 day pulsation period are about 0.3-0.5 mag in 
the V band. The difference between these and the 1 
mag drop in 2019-20 is clear. We do note, however, 
that Betelgeuse has undergone other periods of drastic 
dimming a few times over the last 100 years. Dimming 
events of comparable magnitude are visible in Figure 1, 
for instance, in the mid to late 1980s and arguably in 
the early 1950s. An argument could be made for the 
existence of a 35-40 year dimming cycle, particularly if 
we take into account that the sensitivity of instruments 
has improved considerably in the last few decades. We 
note that this 3-4 decade variation is of the same or- 
der as the suggested rotational period. While this could 
potentially be a manifestation of rotational modulation, 
confirmation will require ongoing observation. 

'The low amplitude and scarcity of adequate compari- 
son stars make visual estimates less sensitive to smaller 
changes from one pulsation cycle to another. Digital 
photometric observations exist for the last three decades, 
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Figure 1. Lightcurve of a Ori assembled from publicly available data compiled by the AAVSO, from 1928 to present, and from 
the SMEI observations. Horizontal axes are marked in both years (UPPER) and JD + 2400000 (LOWER). Grey points are 
visual estimates, blue are V-band photometry, from AAVSO. Red points are the SMEI data. 


but both the quality and quantity varied over time. 
Most of the publicly available data have been archived 
by the AAVSO and provide good coverage from the mid- 
1980s to the early 2000s and from 2010 onward. To fill 
in the gap, we supplement the AAVSO data set with the 
observations taken with the SMEI (Solar Mass Ejection 
Imager) instrument aboard the Coriolis satellite (Jack- 
son et al. 2004). 


2.1.1. SMEI photometry of Betelgeuse 


SMEI was designed to follow Coronal Mass Ejections 
(CMEs) from the Sun, and in order to do this, stellar 
signals must be removed from its images. About 6000 
stellar sources plus the brightest Solar System objects 
were catalogued and then subtracted from the images. 
It was soon realized, however, that the source subtrac- 
tion procedure used by the mission can be processed 
into time series photometry of the brightest stars in the 
Sky, essentially turning SMEI into one of the early space 
photometry missions (Buffington et al. 2007; Hick et al. 
2007). SMEI observed a Ori from early 2003 to late 
2011 with a cadence of 104 mins. Each year, data collec- 
tion was split between the three cameras whose outputs 
needed to be rectified. Yearly systematics arise from 


the changing thermal conditions in each of the cameras 
(Tarrant et al. 2008). Slow degradation of the camera 
sensitivity is also apparent in the data. 

We could not remove the annually repeating instru- 
mental signals directly, as the timescale is on the same 
order as the variation of œ Ori. Therefore, we relied on 
the ensemble photometry of neighboring stars to derive 
common instrumental characteristics. We inspected ten 
nearby bright stars and selected y Ori, € Ori and 32 Ori 
to generate a template for the instrumental signals. We 
calculated a smoothed systematics curve by calculating 
the medians of the combined relative intensity data of 
these three stars in 4-day windows placed around every 
time stamp of o Ori, and for each camera separately. 

The rectified SMEI light curve of o Ori is the result of 
scaling the raw data with the systematics curve and then 
transforming it to magnitudes using mgmer = 10.0 mag 
as the magnitude zero point. However, the passband 
of SMEI is not the V band, therefore requiring that we 
scale and shift the light curve to match the AAVSO data. 
To compute the appropriate scaling, we determined the 
brightness difference for six other stars with M1-2 spec- 
tral class in the SMEI catalog to be my — mswgr © 
0.15 mag. We found that we needed to stretch the am- 
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Figure 2. Detailed plot of the recent photometric data. Blue: 
AAVSO V-band photometry. Red: rectified and scaled SMET 
photometry. 


plitude by a factor of 1.8 to match the V data points. 
We then averaged the raw photometry points into 1- 
day bins. While the shape of the variation could also be 
passband-dependent to some extent, the scarcity of over- 
lapping V data prohibited us from performing a more 
detailed comparison. The final light curve is plotted in 
Fig. 2, along with the AAVSO V-band data. The SMEI 
data confirm that the star did not dim excessively during 
the LSP minimum occurring in 2007-08. 

A sample table of the processed and binned SMEI pho- 
tometry and the scaled V-band values can be found in 
Table 1. Here, we provide formal errors calculated as 
the shot noise from the number of electron counts. 

The photometric light curve reveals a richer set of fea- 
tures than the visual light curve. The SMEI observa- 
tions, in particular, show both the slow variation from 
the LSP along with additional smaller, more rapid varia- 
tions. The SMEI data also put the severity of the recent 
dimming event in perspective: the brightness of the star 
did not drop below 1.1 mag in the V band during the 
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Figure 3. Top: power spectrum of the photometric obser- 
vations. The strongest LSP frequency and the position of a 
yearly alias are indicated. The insert shows the spectral win- 
dow function of the data with the prominent yearly aliases. 
Second plot: spectrum after we removed the LSP signal. 
Pulsation peaks and aliases indicated. Third plot: residual 
spectrum after both the LSP and the pulsation frequency 
removed: fı marks the significant peak that remained. Bot- 
tom plot: power spectrum with the LSP removed from the 
data, in log space: we fitted the granulation noise component 
with a 1/f function and the pulsation frequency region with 
a Lorentzian profile. 


last 40 years, whereas the dip commencing in November 
of 2019 dimmed the star to 1.6 mag in that band. The 
light curve also highlights some smaller variations on the 
order of a few hundredths of a magnitude on timescales 
of days to weeks. Similar variations are present in the 
SMEI light curves of other nearby stars as well, so we 
do not consider these to be an intrinsic feature of œ Ori. 


2.1.2. Frequency analysis of observations 


We analyzed the frequency spectrum of the photomet- 
ric light curve with Period04 (Lenz & Breger 2005). We 
are able to identify the LSP and pulsation frequency re- 
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gions easily, as shown in Fig. 3, but the identification 
of individual frequency components intrinsic to the star 
was hindered by the presence of yearly aliases. Most no- 
tably, the — fusp + 1/yr component coincides with the 
pulsation frequency region. As the FM pulsation mode 
itself is only slightly longer than one year, its harmonics 
and/or overtones could coincide with yearly aliases. 

We first apply a pre-whitening procedure to the data 
with LSP components. Figure 3 shows that the LSP is 
not strictly cyclic and that o Ori hovered in a bright 
state throughout the 2010s. We test combinations 
of multiple harmonics and subharmonics of the main 
fuse = 0.000423 d-! frequency (Prsp = 2365 + 10 d), 
which is 1596 longer than that determined by Kiss et al. 
(2006), but still within their uncertainty range. We 
use the 0.5 and 2.5 frsp components for the final fit, 
which successfully reproduces the deep LSP minima in 
1985/1989 and in 2001/2006-7. Non-sinusoidal features 
in LSP light curves are common for smaller red giants 
in the Magellanic Clouds: one half of the LSP cycle 
shows an eclipse-like dip, and the other half resembles 
a plateau. The model proposed by Soszyński & Udalski 
(2014) to reproduce this shape assumes a nearby orbit- 
ing companion and associated dense cloud. Currently, 
there is no indication of a companion orbiting a Ori, 
but some observations suggest that the recent dimming 
was likely caused by excess dust?. (Levesque & Massey 
2020; Safonov et al. 2020). 

We detect two significant frequency components 
(fpusi = 0.002469, fpuise = 0.002213d ^!) at the pul- 
sation frequency peak, in agreement with the expected 
short lifetime of the mode. We likewise detect the first 
harmonic (2 fpuis) of the stronger pulsation component. 

Since the pulsation signal is non-coherent, we fit it 
with a Lorentzian profile as in Kiss et al. (2006), but in 
combination with a 1/f component to account for the 
red noise component of the convective motions (bottom 
panel of Fig. 3). We calculate a pulsation frequency of 
fouls = 0.00240 + 0.00014 d-! from the peak of the pro- 
file, corresponding to a period of Pous = 416+ 24d. We 
can also use the the full width at half maximum (T) of 
the profile to estimate a mode lifetime of 7 = 1/aT = 
1174d, or & 3 pulsation cycles. The mode lifetime 
matches the value calculated by Kiss et al. (2006); the 
pulsation period, however, is 7.2% longer, though still 
within their uncertainty range. We note that Dupree 
et al. (1987) determined a similar, 420 d period, but 


? We note that this view is not held uniformly, however; e.g. Dhar- 


this was based on only three years of photometric obser- 
vations. 

Apparent changes to the period likely arise from (1) 
the non-coherent nature and short lifetime of the mode 
and (2) interference with photometric variations caused 
by convective motions and the evolution of hot spots. 
Differences of up to 1596 among apparent periods calcu- 
lated from shorter and longer data sets have been found 
for other RSGs as well (Chatys et al. 2019). Presently, 
we report a new period for the photometric data covering 
only the last three decades; disentangling the temporal 
evolution of the pulsation is beyond the scope of this 
work. 


2.1.3. Detection of the first overtone 


After pre-whitening the data using these frequen- 
cies, one significant peak remained at fı = 0.005392 + 
0.000002 d^! (P, = 185.5 + 0.1d)?. Neither this com- 
ponent nor the harmonic was described by Kiss et al. 
(2006), nor is it present in the power spectrum of the 
complete visual light curve. However, f, can be iden- 
tified in some segments. This peak could suggest the 
presence of the first overtone with a period ratio of 
P,/Po = 0.445 + 0.014 (using the Lorentzian fit to Po) 
in o Ori. Overall, we expect the P;/Pp period ratio to 
be ~ 0.5 for red giant and supergiant variables, though 
model predictions typically focus on lower mass ranges 
(Fox & Wood 1982; Kiss et al. 1999). We discuss period 
ratios derived from our own models in Section 4.2. 

Multi-periodicity is not uncommon among red giants 
and supergiants. Kiss et al. (1999) detected more than 
one periodicity in 60% of 93 well-observed, semiregular 
variable stars. Although some of these were a combi- 
nation of pulsation and LSP signals, 30% of the stars 
clearly pulsated in the fundamental mode and first over- 
tone simultaneously. We therefore consider it likely that 
a Ori also pulsates in more than one mode, though we 
cannot conclusively exclude the possibility that the fı 
signal detected in our analysis corresponds to a yearly 
alias or a harmonic of the non-coherent pulsation signal 
that the photometric data do not resolve properly. 

It would be informative to collect photometric obser- 
vations of a Ori throughout the year for as long as pos- 
sible in order to minimize the gaps in the data and di- 
minish such aliasing in the frequency domain. 


2.2. Timing of minima 


mawardena et al. 2020 do not agree) 


3 Uncertainties for fı were calculated with the assumption of a 
single coherent Fourier component: more data will be needed to 
assess the validity of this assumption. 
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A standard means of identifying deviations from an 
assumed periodic signal is the O-C method.* Here we 
attempted to identify and time the various larger and 
smaller minima in the light curve. The light curve data 
appear to alternate between two states: one defined 
by deep minima exceeding 0.5 mag (e.g., at JD 49800, 
52750, 54000, 54400, 58500 and the dip itself at 58800), 
the other by shallower and more frequent meandering 
(e.g., around JD 51500, 53200, 55000 to 57000). How- 
ever, the annual gaps make it difficult to identify enough 
minima accurately, and it is thus possible that we simply 
miss one type during certain intervals. Time spans be- 
tween consecutive shallow minima can be as short as 60— 
100 days—much shorter than the FM pulsation period. 
We see no indication of discrete frequency components 
corresponding to these intervals in the power spectrum 
of the star, which suggests they are not high-degree pul- 
sation modes. The timescales and low amplitudes, how- 
ever, do match the convective turnover times of giant 
convection cells: our photometric results agree with pre- 
dictions of timescales from 3D radiative hydrodynamic 
models and the time-resolved results of spectropolari- 
metric observations of the surface of the star (Freytag 
et al. 2002; López Ariste et al. 2018). We therefore 
attribute these short-timescale, stochastic variations to 
the photometric effects of convective cell turnover on the 
surface of a Ori. 

The critical observational features of Betelgeuse are 
summarized in Table 2. We add a final note that obser- 
vations released since the presentation of the lightcurve 
in this work show a new brightness minimum occurring 
less than half a pulsation period later, but at a regular 
depth of V ~ 1.0 mag (??). These observations are con- 
sistent with behavior seen in our longitudinal data, par- 
ticularly those periods in which the short-lifetime mode 
shifts in phase abruptly. It is likewise consistent with 
the photometric effects of changing convective cells. Our 
results thus demonstrate that these two phenomena are 
sufficient to explain the subsequent, short-term semireg- 
ular variability of a Ori without invoking explanations 
such as multiple, opaque dust clouds orbiting the star. 


3. CLASSICAL EVOLUTIONARY MODELS 


Having carefully collated the set of observational cri- 
teria described in Table 2, we proceed in modeling the 
system. Our numerical efforts include three types of 
simulation: (1) classical evolutionary tracks; (2) linear 
pulsation models; and (3) short-timescale, 1D, implicit 


hydrodynamical evolution. We discuss results from each 
in this order. 

We compute evolutionary tracks for stellar models 
with initial masses of 10-26 Mo,;. Calculations are car- 
ried out from the pre-main sequence to the termination 
of the helium-burning giant branch, with the terminat- 
ing condition set by the amount of helium remaining in 
the core of the star (M(*He) ~ 1078 Mo). Models in 
an evolutionary phase more advanced than core helium 
burning are less favored probabilistically, as the star will 
spend considerably less time in such phases. As shown 
above, they are also unlikely to be consistent with the ex- 
isting array of observational constraints, especially since 
these constraints prove to be discriminating even within 
the set of strictly core helium-burning models. As such, 
we do not consider post-core helium-burning models in 
further detail. 

Figure 4 shows a set of evolutionary tracks evolved 
from the zero-age main sequence (ZAMS) to the end of 
core helium burning. Masses indicated refer to the initial 
mass. In subsequent discussion, we refer to models by 
their initial masses, though typically the mass of the star 
will be between 2-3 Mo smaller at the termination of 
its evolution (and onset of its hydrodynamic evolution) 
due to mass loss during its prior stages. 

Our initial grid of models does not invoke rotation and 
has fixed, solar metallicity represented by a heavy metal 
fraction of Zi, — 0.02. We consider multiple values of 
the convective mixing length owrr, ranging from 1.8 
to 2.5. As massive stars are quite sensitive to the pre- 
scriptions used for convective boundaries and convec- 
tive overshoot, we adopt convective overshoot settings 
of fov, = 0.010H,° surrounding hydrogen- and helium- 
burning zones (Herwig 2000; Paxton et al. 2018). We 
set convective boundaries according to the Schwarzschild 
criterion, and we do not include semi-convection®. We 
use the Cox prescription of the mixing length formalism 
(Cox & Giuli 1968). We do not use the MLT++ op- 
tion in the calculation of our evolutionary models. Our 
surface boundary conditions are set by a simple pho- 
tosphere model, whose implementation is described in 
Paxton et al. (2011). 

We account for mass loss in the evolutionary cal- 
culations via MESA's implementation of the “Dutch” 
wind schemes, a composite of prescriptions summa- 
rized in Reimers (1975); de Jager et al. (1988); Bloecker 
(1995) and van Loon et al. (2005). We model the low- 


5 Multiples of the pressure scale height, dln P/d]n T. 


4 O—C refers to the observed minus calculated method, where we 
measure the time differences between observed events (e.g., cycle 
minima or maxima) and a periodic ephemeris. 


6 only available when the Ledoux criterion (i.e. composition gra- 
dient) is used to set convective boundaries 
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Figure 4. 

(UPPER) A set of classical evolutionary tracks for 10-25 Mo 
computed with MESA. Initial mass per track as indicated. 
All models shown are computed until the end of helium burn- 
ing and shown from ZAMS. All tracks adopt a mixing length 
of amir = 2.1. (LOWER) Same as above, but rescaled to 
highlight the evolutionary region consistent with the tem- 
perature constraints provided by Levesque & Massey (2020), 
shown in pink. The extended temperature regime permitted 
by the mixing length degeneracy is shown in blue. 


temperature mass loss via the prescription of de Jager 
et al. (1988), adopting a wind coefficient of 7 = 0.8 as 
default. 

We test a range of 7 values and find that while the 
choice of 7 does impact the terminal mass of the evolu- 
tionary model, our results are predominantly sensitive 
to the radius. The relationship between an evolutionary 
model’s terminal radius and its input controls—mass, 
metallicity, mixing length, convective overshoot, mass 
loss coefficient, etc.—is complex, and we do not gain 
much insight on this interplay by varying 7. We do 
not use mass loss or rotation during the hydrodynamic 
evolution itself, as the impact of these processes on a 
timescale of several decades is negligible. 

A critical component of our classical modeling objec- 
tive involves reproducing the recently observed temper- 
ature of a Ori. However, given limited a priori informa- 
tion on the star’s mass and evolutionary phase, there is 


a strong degeneracy between choice of amur and pre- 
dicted temperature. While this issue emerges even for 
well-constrained systems (Joyce & Chaboyer 2018a,b), 
the magnitude of the degeneracy is exacerbated as ob- 
servational constraints loosen and the structural com- 
plexity of the stellar models increases. Hence, even if we 
assume that the atmospheric models used by Levesque 
& Massey (2020) can determine the temperature cor- 
responding to the observed line profile with high accu- 
racy and precision, the underlying evolutionary models 
themselves may shift by about +200 K. This introduces, 
at minimum, the same uncertainty on the evolutionary 
stage at which the star crosses the observed tempera- 
ture. It is likewise prudent to extend the uncertainties 
on our temperature measurements anyway, as RSG tem- 
peratures are notoriously difficult to infer. In particular, 
Davies et al. (2013) note that an SED fitting approach 
would give a higher temperature for Betelgeuse than 
the one we adopt here, and these authors raise concerns 
about the validity of the TiO-band fitting method use 
by Levesque & Massey (2020). Likewise, Ireland et al. 
(2008) find that fits based on TiO and made under the 
assumption of LTE can give temperatures that are much 
too low; based on these findings, Betelgeuse could easily 
have an effective temperature that is hotter by 200 K. 

A looser interpretation of the temperature constraints 
is implemented in practice by extending the observa- 
tional boundaries by a theoretical uncertainty of appro- 
priate order. This is done by measuring the shift in 
temperature a track of fixed mass undergoes when its 
mixing length is adjusted to extremal values. In the 
case of our grid, this shift is calculated for aypr = 1.8 
vs &QMLT = 2.5 and corresponds to a shift in mod- 
eled temperature, in the relevant part of the HR dia- 
gram, of roughly 0.1 dex for a track with middling mass 
17 Mo. We likewise note that the mixing length is not 
the only parameter that contributes to uncertainties in 
the derivation of RSG temperatures; others include con- 
vective overshoot and semi-convection, both of which 
also affect the appropriate choice of mixing length itself. 
We do not account for variations in either in this study, 
but refer to Chun et al. (2018) for a detailed analysis of 
the impact of theoretical assumptions on RSG temper- 
atures. 

We choose to account, approximately, for variations 
among temperature estimates caused by differences in 
observational inference methods and theoretical param- 
eter choices as follows: We expand the observational 
temperature constraints by an amount equal to the shift 
in modeled temperature for two otherwise identical stel- 
lar tracks that adopt minimal and maximal values, re- 
spectively, of apr. This adjusted temperature bound- 
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ary is indicated by the blue strip in Figure 4. The ef- 
fective temperature constraints of Levesque & Massey 
(2020) alone are shown in the much more restrictive pink 
band. The observational constraints on luminosity are 
not strong and do not themselves rule out any of the 
models shown in Figure 4. 

Attempts to reproduce Betelgeuse's present-day rota- 
tion of 5 kms"! (vsini = 5.47 + 0.25 kms" !, Uiten- 
broek et al. 1998; Kervella et al. 2018) with single-star 
evolutionary models are unsuccessful. To this end, we 
compute tracks that use an initial surface rotation of up 
to Q = 0.6504, or roughly 200 km/s on the ZAMS, in 
accordance with Wheeler et al. (2017). In cases where 
the models do not fail outright, the results are not con- 
sistent with even the most generous interpretation of the 
observational constraints. Among tracks that converge, 
even those with the highest values for Q; still fail to pre- 
dict a present-day rotation rate in the vicinity of the 
observed value. 

In particular, tracks with initial rotations approach- 
ing breakup velocity (Q/Qerit ~ 0.7) fail to intersect 
the (extended) effective temperature regime with large 
enough present-day surface rotations. The highest val- 
ues attained by our grid only just reach 1 km/s, and 
these correspond to models with initial masses as low as 
6-10 Mo. Such low-mass models are easily ruled out by 
other constraints, especially period. 

Our results from this exercise are thus similar to those 
of Wheeler et al. (2017), who find that “models at the 
tip of the RSB typically rotate at only ~ 0.1 km/s, in- 
dependent of any reasonable choice of initial rotation.” 
Though Wheeler et al. (2017) are able to create rotat- 
ing models consistent with 3e uncertainties on their ob- 
servational constraints at the time, our constraints pri- 
oritize the fundamental mode frequency and include a 
much tighter range on effective temperature. More so- 
phisticated modeling of the rotational aspects of a Ori's 
evolution are beyond the scope of this paper. 

'The terminal models from the evolutionary run pro- 
vide both the structural input for calculations with the 
linear pulsation program (next section) and the initial 
conditions for the hydrodynamic study (Section 5). 


4. SEISMIC MODELS 


Used in conjunction with classical parameters, syn- 
thetic frequencies are an extremely powerful tool for dis- 
criminating among possible models of a star. The case 
of a Ori is no exception. 


4.1. Linear perturbations 


We use GYRE to solve the linearized pulsation equa- 
tions for high-resolution structural models produced 
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Figure 5. Adiabatic p-modes are calculated with GYRE 


for all relevant evolutionary tracks. Periods, in days, of all 
models consistent with the observed temperature constraints 
are shown, coded by color for mass and by marker style for 
mixing length, as indicated. (UPPER) Masses range from 
10-24 Mo at a resolution of 1.0M and mixing lengths range 
from 1.8 to 2.5 at a resolution of 0.1. (LOWER) Masses 
range from 10-24 Mo at a resolution of 1.0, Mo and amur 
is fixed at 2.1. Here, the observed temperature constraints 
adjusted to account for the theoretical uncertainty in onurr. 
All models shown adopt 7 = 0.8 and Z = 0.02. The observed 
seismic constraints from Kiss et al. (2006) are indicated with 
blue horizontal lines. 


during the RSG phase (Townsend & Teitler 2013). The 
GYRE program is based on a Magnus Multiple Shooting 
(MMS) scheme and provides both adiabatic and non- 
adiabatic calculations. We consider only adiabatic re- 
sults in this analysis. Figures 5 through 8 show results 
from these calculations. 

As a track that intersects the observational require- 
ments will typically do so at multiple evolutionary 
timesteps, we can produce several pulsation profiles per 
track. Where the models are compatible, we generate 
synthetic frequency spectra at short intervals. In our 
frequency modeling, we do not restrict to a search for 
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Figure 6. Same as Figure 5, but for the first overtone, O1. 
Uncertainties have been scaled by the O1/FM ratio, yielding 


+13.5 days, shown with red horizontal lines. 


FMs (npg = 1”, | = 0,m = 0) alone; rather, we note 
the modes and values of any frequencies in the vicinity 
of Betelgeuse’s two dominant periodicities. However, 
it is immediately clear from our calculations that peri- 
ods constituting the primary decaying sequences in both 
panels of Figure 5 are fundamental mode periods. 

To account for the theoretical uncertainty on Tog, we 
consider two metrics by which an evolutionary track is 
compatible with the observations. In the upper panel 
of Figure 5, we show the periods of adiabatic p-modes 
versus termination age for a collection of models with 
a range of initial masses and mixing lengths; here, the 
structural models used in the seismic analysis have ef- 
fective temperatures strictly within 3600 + 25 (Levesque 
& Massey 2020). 

In the lower panel of Figure 5, all models use apr = 
2.1, but are checked for consistency against the ex- 
tended, theoretical temperature constraints described in 
Figure 4. In both panels of Figure 5, all masses refer to 


7 The lowest radial order for pressure modes, as defined by GYRE 
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Graph 1: Fundamental mode; variable mixing length 
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Figure 7. Additional parameters, including present-day ra- 
dius and mass, are shown as a function of period. Period 
constraints are shown as green, vertical lines. Models in the 
middle and lower panels of each graph are emphasized if they 
are also compatible with the radial bounds set by the inter- 
section in the corresponding radius vs period graph. (UP- 
PER) A random sample of models spanning initial mass, 
mixing length, mass loss parameter, and degree of helium 
exhaustion, not pre-selected for agreement with any temper- 
ature constraints. (LOWER) All models with owrr = 2.1, 
1) = 0.8, and terminal conditions determined by agreement 
with effective temperature, including theoretical uncertain- 
ties. 


the initial mass of the model and the 416 + 24 day peri- 
odicity is denoted by a blue horizontal band. 

In the upper panel of Figure 5, the period-age se- 
quence is tighter and more well-defined, but the results 
between the upper and lower panels are largely consis- 
tent. We note that sub-sequences comprising stars of 
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Graph 3: First overtone; variable mixing length 
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Figure 8. Same as Figure 7, but derived according to seismic 
agreement with the first overtone, O1. 


particular mass are more apparent in the fixed owrr 
case. This visualization more clearly suggests that, at 
least for masses in the 15-20 Mo range, there will nec- 
essarily be some point along the helium-burning branch 
during which the star will pass through the appropriate 
frequency band. However, the temporal window during 
which this occurs is quite small in the context of evo- 
lutionary timescales-on the order of 0.5-1.0 Myr. The 
requirement that this time frame align with a particular 
observed temperature ends up being quite restrictive. 

Collectively, these results suggest a median, model- 
derived mass range of 16-21 Mo (with some outliers 
as high as 24M), at a resolution of 1 Mọ. This is 
broadly consistent with other modellers’ results, though 
our results are more accommodating at the lower-mass 
end. 


We repeat this analysis using as a period constraint 
the probable overtone observed in our photometry (Sec- 
tion 2.1.2) and confirmed in our synthetic frequency 
spectra: an O1 mode oscillating with a period of 185d. 
In lieu of an independent value, we scale the uncertainty 
of the fundamental mode to derive an uncertainty on 
the overtone of 13.5d, yielding P, = 185 + 13.5. This is 
shown in the red, horizontal bars of Figure 6. Though 
suggested by our photometric analysis, confirmation of 
the detection of this mode and its classification required 
supporting evidence from theoretical spectra. This ar- 
gument is detailed further in Section 4.2. 

Figure 7 shows other fundamental parameters as a 
function of period. The FM and its uncertainties are 
defined by green, vertical bars in all panels. An ana- 
log using the O1 period and its uncertainties, defined by 
dark red vertical lines, is shown in Figure 8. 

Models in the upper panels of Figures 5 and 6 span the 
full range of masses and mixing lengths considered in our 
grid and additionally vary in the prescribed mass loss co- 
efficient (7) = 0.2-1.0), but they are not pre-restricted by 
agreement with temperature constraints. Instead, these 
evolutionary tracks are terminated at arbitrary intervals 
along the helium-burning branch, with spacing set by 
the degree of helium exhaustion in the core. This is done 
to produce a more well-populated sequence that incorpo- 
rates additional sources of uncertainty in the modeling 
assumptions. Despite this added theoretical noise, the 
range of possible radii across all models remains heavily 
restricted by the observational period constraints. 

All models in the lower panels of Figures 5 and 6 inter- 
sect the theoretically extended temperature uncertain- 
ties (which essentially sets their termination ages) and 
adopt ary = 2.1 and 7 = 0.8. 

In the uppermost panel of Graph 1 in Figure 7 and 
Graph 3 in Figure 8, we show radius as a function of 
FM and O1 period, respectively. Though there is some 
scatter in the synthetic data, the radial span of the 
period-compatible models is very narrow in both cases, 
especially compared to the range of radial estimates col- 
lated in Dolan et al. (2016). We recall from earlier dis- 
cussion that literature estimates of Betelgeuse's radius 
range from ~ 500R to nearly 1300 Ro, whereas the 
results presented here suggest little possibility outside 
700 — 900 Ro. If we interpret the FM period mea- 
surements as hard limits, our results suggest a radius 
for a Ori of 764 Re, with lo uncertainties of roughly 
30 Ro and non-symmetric limits of Rmax = 880 Re and 
Rmin = 702 Ro. Figure 8 suggests consistent but even 
tighter limits of approximately 700 — 800R;; We thus 


report a 3c, model-derived radius of 764." Ro. 
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mode periods computed with GYRE is shown for the same 
models and temperature definition used in the upper panel 
of Fig. 5 (variable aurr). The range indicated by horizontal, 
green lines is P1/ Po = 0.44 + 0.04. 


In the middle and lower panels of each graph, we show 
the models’ initial masses and terminal masses, respec- 
tively, as a function of period. In these plots, we em- 
phasize those models that also have radii within our 3c 
uncertainty bounds with larger, darker markers. Con- 
sidering all possible observational constraints, we report 
model-derived estimates for the initial and present-day 
masses of Betelgeuse as approximately 18-21 Mo and 
16.5-19 Mo, respectively. Though best mass ranges dif- 
fer slightly among them, these represent median values 
collated from the four variants of this test (Figures 7 
and 8), encompassing two interpretations of the mix- 
ing length-temperature degeneracy for each of Po, P. 
The lower rows of Table 2 summarize these findings. 
Taking into account the likelihood of a previous merger 
event, which would significantly complicate any infer- 
ences about the state of Betelgeuse at birth, it is our 
present-day mass estimates that are most pertinent. 


4.2. The first overtone 


In Sect. 2.1.2 we presented observational evidence of a 
new frequency component, fi, that corresponds to a pe- 
riodicity of P, — 185 d. While strong aliasing caused by 
annual gaps in the data makes this detection somewhat 
uncertain, the ratio of this mode to the fundamental sug- 
gests that it is the first overtone. Because the literature 
on seismic models of RSGs is sparse (with most works 
focusing instead on the low- to medium mass regime; 
c.f. Trabucchi et al. 2019), we supplement the photo- 
metric analysis with theoretical results from GYRE. To 
test whether the seismic models agree with the observed 
period ratio we inspect GYRE's prediction for the fre- 
quency of the P, (ny = 2, | = 0, m = 0) mode. With the 


uncertainties propagated from each period, the relevant 
range is P/P = 0.445 + 0.041. 

We find that GYRE’s prediction for Pj is strongly con- 
sistent with 185 d and that the observed period ratio 
gives mass and age predictions that are self-consistent 
with the mass and age ranges constrained by the funda- 
mental mode, as demonstrated in Figure 9. When mix- 
ing length is varied and temperature is fixed, the age 
and initial mass ranges inferred from the period ratio 
are 7-11 Myr and 16-23 Mo, respectively. 

If we consider the progression of the models, we find 
that the period ratio drops in even younger, higher-mass 
models, with P/P as small as 0.30-0.36 at around 
5 Myr. Older, lower-mass models reach the 0.50-0.55 
range, in agreement with period ratios observed in 
semiregular and Mira variables (Kiss et al. 1999; Molnár 
et al. 2019). We thus conclude that the signal detected 
in our photometric data is indeed consistent with the 
first overtone, making Betelgeuse a double-mode star. 

Observing multiple modes in a pulsating star can pro- 
vide stringent constraints on its physical parameters: 
this is the main principle behind asteroseismology (Aerts 
2019). The Fourier analysis in Sect. 2.1.2 only provided 
us with a formal error for the single detected frequency 
peak, but we expect the overtone to have a limited mode 
lifetime, just like the fundamental mode. As Figures 6 
and 8 show, and as is discussed in Section 4.1, the first 
overtone prefers largely the same models as the funda- 
mental mode. 


4.3. Seismic parallax and luminosity 


With the radius of o Ori heavily constrained by the 
seismic models, we can calculate the distance to the 
star based on the measured angular diameter. Us- 
ing an angular diameter of 42.28 + 0.43 mas and the 
30 uncertainty range of the seismic radius estimate, 
764* 116 Ro, we calculate 1687! pc for the distance and 
m = 6.06*025 mas for the parallax. Our values are in 
agreement with the parallaxes derived entirely or in large 
part from the Hipparcos measurements (see van Leeuwen 
2007 and Harper et al. 2008), and place a Ori nearer to 
us. It is, however, somewhat in tension with the more 
recent results based on radio observations, with disagree- 
ment at the 1—2c level (Harper et al. 2017). Figure 10 
shows our results in context. 

This discrepancy could stem from various observa- 
tional or theoretical shortcomings. One possibility is 
that the period shift is caused by large-amplitude, non- 
linear pulsation. Stellar structure adjusts dynamically 
to the changes caused by coherent pulsation, which may 
cause a shift in the eigenfrequencies of the structure. 
Therefore even if the physical parameters of a linear seis- 
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Figure 10. A summary of recent literature distances to o Ori. 
The gray bar represents the seismically-derived values re- 
ported in this work. 


mic model agree with those of the star, the calculated 
and observed periods may not. In the case of p-modes, 
the radius relates to the pulsation period as R ~ P?/? 
for a given mass. From this alone, we estimate that the 
linear period of Betelgeuse should be 500 + 40 d if its 
radius is 887 Ro, as adopted by Dolan et al. (2016). 
This means that a ~ 20% non-linear decrease would 
be required to reproduce the observed 416 + 24 d pul- 
sation period. Given that we cannot yet evolve RSG 
hydrodynamic models to full-amplitude pulsation, we 
cannot infer the period shift directly. For comparison, 
studies show that pulsating Mira models produce pe- 
riod shifts of up to -23% and +15%, which is of the 
appropriate order (see, e.g., Lebzelter & Wood 2005; 
Ireland et al. 2011). However, non-linear period shifts 
scale with the pulsation amplitude. As such, the rela- 
tively low-amplitude, short-lifetime mode seen in œ Ori 
makes such a large shift implausible. Thus, non-linear 
pulsation could be, at best, only partially responsible for 
the discrepancy we observe. 

Another possibility is that the true Rosseland angular 
diameter of the star is smaller than the 41.8 mas value 
adopted in the present analysis, and thus the diameter of 
the photosphere could be as small as 36 mas. However, 
this would suggest that none of the direct imaging and 
interferometric observations, including the multi-band 
models created by Perrin et al. (2004) and Montargés 
et al. (2014), were capturing the photosphere itself, but, 
rather, only two separate layers in the extended atmo- 
sphere. We consider this very unlikely. 

On the other hand, consistency between our parallax 
estimate and the Hipparcos value could indicate that 
the astrometric cosmic noise has been underestimated in 
previous studies, or that it is correlated over the times- 
pan of the astrometric observations. 


Harper et al. (2017) illustrates that the addition of 
a large cosmic noise term to the radio data would be 
required to bring the combined Hipparcos+radio par- 
allax value close to 6.0 mas. Nevertheless, this could 
be an indication that various effects, such as photo- 
center displacements and non-axisymmetric stellar disk 
shapes caused by convective hot spots, have been previ- 
ously underestimated. A resolved image obtained from 
ALMA shows a large hot spot towards the disk limb 
with a temperature contrast of AT ~ 1000 K, which co- 
incides with the rotational axis but provides no context 
for the temporal evolution of the photocenter displace- 
ments (O’Gorman et al. 2017). Kervella et al. (2018) 
proposed that this hot spot corresponds to a rogue con- 
vection cell that might also be magnetically connected to 
ongoing mass ejection from the polar region of the star. 
Presence of one or more persistent spots could mean 
that the cosmic noise level of the star is higher than re- 
ported by Harper et al. (2017), and/or that it needs to be 
modeled as a correlated noise source. Confirming or rul- 
ing out this possibility would require either a sustained, 
multi-year interferometric observational campaign or the 
development of accurate 3D physical simulations of su- 
pergiant stellar atmospheres that can handle the evolu- 
tion of hot spots. Finding new ways for Gaia to process 
heavily saturated stars and thus allowing us to obtain 
better parallaxes and possibly astrometry of individual 
hot spots is yet another avenue forward (Sahlmann et al. 
2018). However, any of these endeavors would require 
substantial investment from the respective experts. 

Finally, with a seismic parallax in hand, we can esti- 
mate a tighter luminosity range based on the Stefan- 
Boltzmann law. Adopting the strict T.g uncertainty 
range reported by Levesque & Massey (2020), we re- 
port the luminosity constraints of a Ori to be logio L = 
4.94* 006, where L is in units of Lo. The effective tem- 
perature range permitted by taking into account the- 
oretical uncertainties extends this range to logio L = 
4:940 0. We note that if we were to superimpose this 
range on Figure 4, it would intersect the RSB at the 
appropriate temperatures for tracks with initial masses 
between 16-21 Mo. This range is self-consistent with 
our other means of estimating mass; however, as mass 
is a derived quantity, we do not use it to restrict the 
domain of our seismic grid. 


5. HYDRODYNAMIC ANALYSIS 


The third component of our modeling relies on 
MESA’s implicit hydrodynamics solver, which we use to 
explore the non-linear oscillatory behavior of the mod- 
els’ envelopes on decadal timescales. Though the term 
“non-linear” can accurately be used to describe any use 
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of the Euler equation, which is non-linear by nature, or 
to describe any models that make full use of the stellar 
structure and evolution equations as opposed to pertur- 
bation analysis, we apply this term only to particular hy- 
drodynamic behavior. We use the term “linear” in the 
hydrodynamic context to describe any oscillation that 
does not excite other modes or cause the development 
of shocks. 


5.1. Method 


We use MESA (version 8118—Paxton et al. 2015) to 
solve the stellar structure and evolution equations by 
means of an implicit hydrodynamics scheme in the La- 
grangian formalism. Critically, this means that velocity 
becomes a dynamical variable similar to density, lumi- 
nosity, and other physical quantities. 

To capture the shock propagation, artificial viscosity 
q is included in MESA’s hydrodynamical scheme (Richt- 
myer & Morton 1967). However, in our calculations, the 
global motion of the H-envelope remains subsonic (~ 0.5 
km s^!) until the end of the simulation. Since this is 
lower than the typical speed of sound in the atmosphere 
(~ 10 km s^1), capturing shocks is not integral to our 
work, and hence q acts only as a safeguard to prevent 
numerical instabilities. 

Furthermore, the code uses the energy-conserving, 
time-discretization scheme of Grott et al. (2005), which 
ensures that the models at two consecutive timesteps are 
consistent with each other. We describe the precise con- 
figuration for our simulations in more detail in Appendix 
A, and provide the inlists necessary for reproduction?. 

It is important to note that MESA’s implicit hydrody- 
namics scheme does not include the equations of stellar 
pulsation directly; rather, MESA has a dedicated mod- 
ule, RSP (Radial Stellar Pulsations; Smolec 2016), for 
calculating non-linear pulsating models. However, RSP 
is not suitable for stars with the luminosity-to-mass ra- 
tios (L/M) of giants such as a Ori (Paxton et al. 2019). 
Previous work has shown that the pulsation modes of 
stars with L/M ratios similar to that of Betelgeuse have 
high enough growth rates to induce shocks even if the 
implicit solver is employed (Heger et al. 1997; Paxton 
et al. 2013; Smolec 2016; Yoon & Cantiello 2010; Gold- 
berg et al. 2020). 


5.2. Period-Radius estimates from hydrodynamic runs 


While the main goal of incorporating hydrodynamic 
simulations into our analysis is to study a canonical 
model of Betelgeuse's envelope rigorously, non-tailored 


8 We will publicly release the inlist files on Zenodo at the conclusion 
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Figure 11. We show the short-timescale, hydrodynamic 
evolution for two MESA models with observationally consis- 
tent features at the termination of their evolution in terms 
of normalized luminosity (UPPER) and normalized radius 
(LOWER) versus time. From these simulations, we can ex- 
tract cycle lengths as a secondary means of estimating oscil- 
lation periods. 


hydrodynamic simulations also provide a second method 
of calculating short-order pulsation modes, and thus a 
means of independently verifying the linear calculations. 
From a small grid of cursory hydrodynamic runs of vary- 
ing mass, we can estimate theoretical pulsation cycle 
lengths directly. 

We first conduct an exploratory investigation of the 
hydrodynamic evolution for a small subset of the models 
in our grid, restricting to those with initial masses be- 
tween 17-23 Mo. We require that the timestep does not 
exceed some fixed, small value—typically 5000-10,000 
seconds—and compute the temporal evolution for sev- 
eral decades. 

Figure 11 demonstrates the oscillatory behavior of 
two hydrodynamic models with slightly different ini- 
tial masses. If not handled correctly, the hydrodynamic 
models will rapidly expand from their evolutionary ini- 
tial conditions—in most cases, to nearly double our ra- 
dial limits—before stable pulsations emerge. This is 
caused by a discrepancy between the luminosity of the 
inner boundary of the simulated stellar envelope and 
the actual stellar luminosity. Thus, over time, the star 
deposits its energy near the surface, making the star 
expand. This can be mitigated by applying relaxation 
procedures to the initial hydrostatic model (Wood et al. 
2004; Nicholls et al. 2009b; Ireland et al. 2011; Saio et al. 
2015). 

However, it is still possible to derive the pulsation pe- 
riods and average radii of these models based upon se- 
lected cycles before shocks and/or numerical failure oc- 
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Figure 12. The left panel shows the sequence of stellar radii against pulsation periods, extending from RR Lyrae and Cepheid stars 
to Miras and RSGs. The right panel highlights the region containing observations of other variable red giants and measurements 
extracted from the hydrodynamic simulations, demonstrating that Betelgeuse's 416 day, rather than 2050 day, periodicity lines 
up better with the modeled sequence. In both panels, gray dots correspond to observations of lower-mass pulsators. Variable 
red giants in particular are shown in open red circles, with Betelgeuse's two modes represented by closed red circles. Points 
derived from models are colored, with the colorbar indicating their mass. 


cur, thus providing a cursory but independent validation 
of the pulsation periods computed with GYRE. 

'The modeled data are produced by estimating the in- 
stantaneous period and radius values from the hydrody- 
namical models using a combination of quadratic and 
a sine functions fit to short segments of the radial evo- 
lution. In this way, we can extract multiple theoretical 
R, P estimates from one hydrodynamic model. We com- 
pare these to a set of direct and inferred observations 
of pulsation periods and radii of variable stars. The 
bulk of these data come from the collection of Szatmáry 
(2004),° which contains a variety of variable stars in- 
cluding RR Lyrae, Cepheids, and Miras. In addition, 
we collate period and radius estimates for a number 
of other supergiants from available literature: CE Tau, 
TV Gem, o Sco, V766 Gem, AH Sco, UY Sct, V602 Car, 
VY CMa and KW Sgr (Wasatonic & Guinan 1998; 
Levesque et al. 2005; Kiss et al. 2006; Ohnaka et al. 
2013; Wasatonic et al. 2015; Wittkowski et al. 2017). 
Of these, pulsation and LSP period measurements are 
available for a Ori and TV Gem, and KW Sgr has two 
distinct radii published—these have dashed lines con- 
necting their measurements. 

Figure 12 shows a set of synthetic periods and radii de- 
rived from the hydrodynamic grid. Also shown in Figure 
12 are (1) the P-R sequence constructed from observa- 
tions of pulsators across a wide mass range (gray dots); 
(2) the observed FM and LSP periodicities for the small 


9 http:/ /astro.u-szeged.hu/P-R. relation/pr. poster.html 


number of red giants listed above (red, open circles); 
and (3) the 416 and 2050 day periodicities of Betelgeuse 
(red, closed circles). Masses of the synthetic stars are 
indicated via the color bar. 

It is well-known that acoustic modes scale with the 
average density of the star, which itself largely depends 
on the radius. We should therefore expect a clear corre- 
lation between radius and period, as seen both here and 
in the linear seismic analysis (Figure 7). 

As is clear in the left panel of Figure 12, there is a 
well-defined P, R sequence spanning RR Lyrae up to 
synthetic supergiants. The periods and radii extracted 
from the hydrodynamic models extend the established 
sequence of pulsating stars to higher radii in a system- 
atic and continuous way, indicating that our models ex- 
perience p-mode pulsation until our simulations are ter- 
minated or the numerics break down. A track of given 
mass can produce multiple data points by sampling at 
different times, and hence different degrees of radial ex- 
pansion, during the hydrodynamic calculation. 

In turn, the right panel hints at certain mode clas- 
sifications for some of the observations. The periods 
for a number of stars with measured radii fall cleanly 
on the model sequence, and some fall above: the latter 
could suggest either pulsations in an overtone or that 
their radii have been overestimated. Given the compli- 
cated circumstellar environment surrounding many su- 
pergiants, and our own findings on the radius of a Ori, 
the latter is a plausible explanation. Finally, the LSP 
signals are clearly separate from the model sequence, 
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confirming once again that the 2050 day periodicity is 
not driven by acoustic variations. 

Even before more careful modeling of the hydrody- 
namic evolution, it is evident that pulsation periods 
emerging naturally in the simulations are of the same 
order as Betelgeuse’s 416 day periodicity. The linear 
and hydrodynamic seismic analyses both demonstrate 
that this is œ Ori's fundamental mode, a fact which, 
when combined with other classical observations, places 
particularly strong constraints on the star’s radius. 


5.3. Possibility of self-excitation due to non-linear 
effects 


As stars evolve across the HRD, they may undergo 
mode transitions when a new pulsation mode becomes 
unstable. At this point, the star can switch to the new 
mode—a phenomenon observed directly in RR Lyrae 
stars—or transition to a multimode pulsator. 

Mode growth rates have various definitions. In the 
linear framework, they usually represent the natural 
timescale of changes in the pulsation energy of the star 
(Catelan & Smith 2015). Growth rates are sometimes 
calculated directly from the change in amplitude be- 
tween successive cycles, but one must keep in mind that 
in non-linear calculations, amplitudes do not grow indef- 
initely. Hence, non-linear growth rates only agree with 
the linear values initially, eventually fading back to zero 
(Yoon & Cantiello 2010). Normalized growth rates are 
thus scaled with the pulsation frequencies. In the case 
of, e.g., OGLE-BLG-RRLYR-12245, this mode transi- 
tion lasted for hundreds of cycles, as is consistent with 
the small growth rates of the modes (Soszyński et al. 
2014). But, in contrast to classical pulsators, semireg- 
ular stars have very large growth rates that can lead 
to strong mode interactions, some of which may even 
become chaotic (Buchler et al. 1996). As such, it is 
theoretically possible that Betelgeuse has recently ex- 
perienced a rapid mode transition, or a rapid increase 
in amplitude of an overtone mode already present, and 
that the superposition of the resulting modes created 
the unusually low brightness minimum seen in Novem- 
ber of 2019. It is thus worth investigating whether such 
a situation can be simulated; however, modeling multi- 
mode pulsation in the non-linear regime is notoriously 
difficult; at present, this is only reliably reproducible for 
stars with much lower L/M ratios (Kolláth et al. 2002; 
Smolec 2016). It is thus beyond the scope of the cur- 
rent paper to investigate such a situation, though this 
scenario is one we hope to address in a subsequent in- 
vestigation. 


5.4. Analysis of Canonical Hydrodynamic Model 


We consider the evolution of a canonical hydrody- 
namic model whose initial and terminal evolutionary 
conditions are consistent (as best as possible) with the 
parameters reported in Section 4. We construct a star 
with initial mass 21 Mọ and terminal mass of 19.54 
Mo during core He-burning. In order to force the stabi- 
lized radius to be consistent with our reported values, we 
must inflate the mixing length parameter to amr = 3.0. 
However, we still require that our hydrodynamic model 
intersect the theoretical temperature uncertainties de- 
scribed above, 3600 + 200 K, during its oscillations. 

Following the general methodology outlined in Gold- 
berg et al. (2020), we cease tracing the evolution of the 
innermost 6.0Mo, representing the core, at the conclu- 
sion of the classical evolutionary run. At this point (of 
ten called “core removal” or, more accurately, “core freez- 
ing”), a typical timestep in the code is still 100 - 1000 
years. As a numerical test, we run an evolutionary (non- 
hydrodynamic) model of a 20Mo star and find that after 
He exhaustion, it takes approximately 50,000 years for 
the luminosity to change by 0.05, whereas our hydro- 
dynamic models show luminosity variations of similar 
scale over the course of 40 years. Hence, the core can 
be considered unchanged to good approximation over 
the duration of the hydrodynamical simulations. It is 
thus valid to adopt constant inner boundary conditions 
set by reut such that M(reut) = 6.0Mo, L = L(reut) 
throughout our calculations of the short-timescale, dy- 
namical evolution of the envelope. The value of 6.0Mo 
is chosen so that 1.0 Mo of the outer He-layer remains 
along with the entire H-envelope. The He layer, which 
sets the base of the hydrodynamical simulation, forms 
a core-envelope structure with the H-envelope, and the 
higher density core ensures that the oscillation of the 
envelope does not interact directly with the mass gap. 

In order to maintain a stable configuration after ceas- 
ing to evolve the core, we allow the model to settle into a 
hydrostatic approximation before turning on the hydro- 
dynamic solver. To capture the short timescale motion 
to adequate resolution, we limit the timestep of the hy- 
drodynamical evolution to a maximum of 104 s. A larger 
timestep of only ~ 10° s can already result in the era- 
sure of modes with a sub-annual period; this is due to 
the implicit nature of the hydrodynamical solver. We 
note that an implicit hydrodynamic scheme is necessary 
in order to follow the global motion of the star because 
the relative distances among mass shells near the surface 
are small. In terms of the Courant timescale, it is ~ 108 
larger than that required by explicit time discretization 
(~ 0.5 s). Thus, in order to track the motion of the 
surface with sufficient temporal resolution, implicit hy- 
drodynamics must be employed. 
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Our canonical model is evolved for a total of ~ 10,000 
steps from the initiation of hydrodynamics until the 
point at which stellar expansion begins to disturb the 
pulsation frequency. We find that beyond ~ 30 years, 
the motion in the star becomes large enough to interact 
with the convective layer, causing the timestep to drop 
as low as 1 x 10? s before the code is unable to evolve 
the model forward in time. At this point, we stop the 
simulation. 


5.4.1. Global Features 


We first discuss the temporal evolution of the critical 
observables in the canonical model. In the four panels of 
Figure 13, we plot the luminosity, effective temperature, 
radius and surface velocity of the star. 

'The system enters into a state of pulsation with steady 
but growing cycles a few years after the hydrodynamic 
solver is switched on. Early in the evolution, quasi- 
annual oscillatory behavior is present in the luminosity 
and effective temperature, and the stellar radius exhibits 
a consistent periodic motion on top of a steady expo- 
nential expansion. In this work, we will consider the 
stellar pulsation only when the motion remains linear; 
we note that once the behavior becomes non-linear, the 
timestep becomes too small to follow the pulsation effec- 
tively. Moreover, non-linear pulsations greatly disturb 
the profile of the stellar envelope, particularly in terms 
of opacity and free electron fraction, which makes direct 
comparison difficult. 

In Figure 13, four vertical, dotted lines indicate mo- 
ments at which we compare the instantaneous values of 
the four quantities. The red and green lines correspond 
to timesteps where the luminosity is at a local maximum 
and minimum, respectively. The blue and purple lines 
correspond to timesteps where the surface velocity is at 
a local minimum and maximum, respectively. 

When the star reaches its brightest point in the pul- 
sation cycle, the effective temperature also reaches its 
maximum. Concurrently, the star is in its most con- 
tracted state (radial minimum) and displays a nearly 
zero surface velocity. This is consistent with the be- 
havior of a classical harmonic oscillator where the dis- 
placement is largest during one cycle. Conversely, the 
luminosity and effective temperature are minimal when 
the star is most radially extended, and when the star 
is maximal in surface velocity, the luminosity, effective 
temperature, and radius are near their average values. 
We then observe that as the star continues to expand, 
a wobble in its motion emerges. As indicated in the hy- 
drodynamic evolutionary tracks of Figure 11, the star 
will eventually approach a state in which it is capable 
of developing a hydrodynamical instability in the outer 


envelope—the exact pulsation properties are themselves 
a function of evolutionary stage and mass, as shown 
in Yoon & Cantiello (2010). The evolutionary time at 
which we launch our hydrodynamic simulations is chosen 
so as best to reproduce the constraints on the classical 
parameters found in the previous analysis. We adopt 
a uniform starting point at which the He mass fraction 
in the core is 107^. Though this is somewhat late in 
the helium burning phase and thus does not correspond 
to the evolutionary phase statistically preferred by our 
classical models, this starting condition lends itself to 
more stable hydrodynamic models. We find that when 
the simulations are initiated at earlier or later burning 
phases, it is difficult to produce stable, suitable stellar 
models that fit the radius and luminosity constraints 
derived in previous sections. As our main priority is 
to match these parameters, we do not further investi- 
gate the effects of evolutionary starting condition. We 
require only that our hydrodynamical initial conditions 
yield models that reproduce the radius, luminosity, and 
mass of the evolutionarily preferred models, as these are 
the features relevant for studying pressure-driven pulsa- 
tion in the envelope. 

'The expansion of the outer radius gradually affects the 
P-R relation, as the sound speed travel time increases 
with increasing distance. Analysis of the radius is addi- 
tionally complicated by (1) how the outermost bound- 
ary of the star is defined and (2) radiation pressure out- 
side the photosphere. A rigorous treatment of radiative 
transport is necessary in the photosphere regime, and so 
we stop the simulation to analyze the motion only when 
the radius is beneath this threshold. 

The effective temperature 7;g and luminosity L be- 
have similarly to radius in the simulations. In the for- 
mer case, this is because MESA calculates the effective 
temperature directly from the luminosity and radius via 
Tå; = (L/Ano pg R?), where og is the Stefan-Boltzmann 
constant. Given the slow change of the radius, Teg pri- 
marily mirrors the fast-evolving L. 

Regarding the evolution of luminosity, we note that 
from year 15 onward, the star exhibits periodic motion 
in its brightness. The early motion is highly regular: 
as in the preliminary grid of hydrodynamic models (see 
Figure 11), we observe a quasi-annual rise and fall— 
the correct timescale for the FM. Near the end of the 
simulation, the large oscillation begins to trigger non- 
oscillatory motion in all quantities, and this is responsi- 
ble for the rapid drop in luminosity. 

In Figure 14, we plot the structural profiles of six 
quantities at points indicated by black circles in Fig- 
ure 13. The global features of the density and temper- 
ature profiles show that outermost 1096 of the stellar 
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time (year) 


Figure 13. (UPPER) The temporal evolution of the luminosity (L), effective temperature (Teg), radius (R) and surface velocity 
(Vsurt) for the characteristic model. The red, green, blue and purple vertical dashed lines are added for comparing the four 
quantities at the same time slices. The black dots are the time moments where the stellar profiles are plotted in Figure 14. 


mass has a relatively low density, sitting between 107°- 
10-7 g cm~’, while the temperature lies between 10?5— 
10° K. A small density bump appears at log;o9(1 — q) = 
log;9(1 — M(r)/M) = —3, which is accompanied by a 
sharp fall in temperature. This occurs in order to main- 
tain hydrostatic equilibrium. We note that a density 
inversion can occur only when convective mixing is in- 
efficient. 

It should be noted that a density inversion, in gen- 
eral, cannot exist when there is convection: the mixing 
will smooth out the density difference. The efficiency 
of this process depends on the ratio of the timescale of 
photon diffusion to the dynamical time, as discussed in 
Jiang et al. (2015). A slow photon diffusion time, as 
applicable to the stellar interior, implies efficient con- 
vection. Hence the density discontinuity is an evanes- 
cent phenomenon. However, near the surface where the 


photon diffusion timescale is much shorter—similar to 
what we have observed—inefficient convective mixing 
cannot remove the density inversion. As described in 
Jiang et al. (2015), this density inversion can form and 
be destroyed repeatedly in the optically thin regions. 
Numerical schemes, such as the inclusion of a photon 
“porosity” by modifying the photon acceleration terms, 
can suppress this phenomenon in optically thick regions 
(Paczyński 1969). 

When convection is less efficient and has a timescale 
comparable with the dynamical timescale, the density 
and pressure gradients can change signs, creating a 
Rayleigh-Taylor instability. Through mixing, the ex- 
cess density can gradually reduce via diffusion. How- 
ever, modeling this phenomenon would require a de- 
tailed, time-dependent convective scheme, which is not 
included in this work. For our case, the density in- 
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Figure 14. Analysis of model star's structure at three points selected during the pulsation, indicated by black markers in Figure 
13. LEFT: The density (top panel), temperature (middle panel) and velocity (lower panel) profiles for the moments at the 
luminosity minimum (black solid line), midpoint (red dotted line), and maximum (green dashed line). RIGHT: Same as the left 
panel, but for the luminosity, opacity and free e^ fraction. In both panels, m represents the enclosed mass, M (r). 


version plays a less important role in the luminosity 
evolution, given that this quantity remains steady from 
log;g9(1 — q) = —1 upward (see subsequent discussion). 

The free electron fraction shows that up to log49(1 — 
q) = —3, the matter is partially ionized. Beyond that, 
the low temperature causes the nuclei to recombine with 
the free electrons. The opacity profile is richer; rather 
than falling monotonically like the free electron fraction, 
we see two major opacity bumps near log,,(1—q) = —1.2 
and —2.5. These correspond to the partial ionization 
zones of H-Hel and Hell, respectively (Cox et al. 1973; 
Kiriakidis et al. 1992). 

The velocity and the luminosity vary dynamically dur- 
ing the pulsation. When the stellar luminosity is mid- 
phase, the whole envelope is contracting with a constant 
velocity of ~ 0.7 km s71. The whole star contracts more 
slowly when it is close to its luminosity maximum or 
minimum. Meanwhile, the luminosity profiles show that 
when the star is at its local minimum, the luminosity 


near the interior of the envelope is lower. The oppo- 
site applies during its local maximum. We note further 
that the motion vanishes approaching the core-envelope 
interface. This suggests that the observed pulsation is 
a collective motion of the envelope without interaction 
with the core. 

'To further quantify changes in the stellar profile dur- 
ing pulsation, we plot in Figure 15 the ratio of the adi- 
abatic temperature gradients for the profiles at lumi- 
nosity minimum and medium with respect to luminos- 
ity maximum. In Figure 15, we observe a sinusoidal- 
like variation of the profile with about 7 nodes, located 
at logigg = —1,—1.5, —1.9, —2.2, —2.3, —2.9 and —3.5. 
'The star at luminosity maximum shows the highest tem- 
perature gradient near the surface at logo q ~ —2.2 and 
—3.1. These correspond to the minima in the temper- 
ature ratio profile and the extrema in the opacity ratio 
profiles. 
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Figure 15. Similar to Figure 14, but for the ratio of the 
adiabatic temperature gradient for the profiles at luminosity 
minimum (blue solid line) and luminosity median (purple 
dashed line), with respect to that at luminosity maximum. 
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Figure 16. Kinetic energy as a function of time for the hy- 
drodynamical model. The blue line is an exponential fit of 
the form 1.2 x 104° exp(0.2t), with time in years. 


These trends are indicative of the &-mechanism, and 
thus explain why the pulsation gradually grows over 
many periods. In particular, during contraction, the 
higher opacity prevents the heat from being stored in 
the deeper layers, which in turn prevents unstable en- 
ergy extraction by ionization. 

From this collection of profiles, we can deduce that the 
~ 1 year variation is driven by the collective expansion 
and contraction of the recombined hydrogen layer. We 
thus conclude that it is this &-related interaction driving 
the fundamental pulsation mode in Betelgeuse. 


5.4.2. Literature Comparison 


To examine the growth of the oscillation further, we 
plot the total kinetic energy of the system against time 
in Figure 16. The kinetic energy, which is dominated 


mostly by the motion of the atmosphere, is much smaller 
than both the total energy and the total gravitational 
energy, which are, on average, —6 x 104° and —1 x 10°! 
ergs, respectively. These are several orders of magnitude 
larger than the kinetic energy, which is 1072? erg, in 
agreement with earlier results (see Cox 1980). We note 
that it is the outermost q — 0.1 — 1 that contributes to 
the atmospheric behavior, including the pulsation. This 
corresponds to ~ 2x 104” erg when this matter is moving 
at a speed comparable to the escape velocity. 

We provide an exponential fit in blue on Figure 16 
to characterize the rate of growth of the kinetic en- 
ergy. A function of Ekin = Aexp(btyear) with param- 
eters A = 1.2 x 10? erg and b = 0.2 provide a good 
fit to the hydrodynamic component of the simulation. 
Naively, this suggests that when the oscillation begins 
to grow, we should expect that the outermost layers of 
material will be expelled by the pulsations within a time 
of ~ 83 years. In reality, however, we see this rate level 
off—see, for example, Cox et al. (1966), who showed 
this in some early pulsation models. This is because vis- 
cous and turbulent dissipation limit the maximum am- 
plitude of the star in the non-linear regime. Historically, 
the level of dissipation in 1D pulsation models has been 
tuned to match the observed amplitudes of RR Lyrae 
and Cepheid stars and to reproduce double-mode pulsa- 
tion in the models. This dissipation was first applied via 
the “artificial viscosity” term and later through the eddy 
viscosity and other o parameters of time-dependent, tur- 
bulent convection (see, e.g., Buchler 1990; Takeuti et al. 
1998; Kolláth et al. 1998; Smolec & Moskalik 2008). 

We note that the models in Fig. 17 stop at about +2% 
luminosity variation or less, whereas Betelgeuse itself 
varies by +10-30% in V and about +20-30% in near-IR, 
(the latter being more closely representative of Betel- 
geuse's bolometric variation). As such, there is plenty 
time remaining for the kinetic energy and amplitude to 
grow and eventually saturate at that level, but this is be- 
yond what can be achieved with hydrodynamics today 
before encountering a numerical runaway episode. 

We note that when the oscillation becomes strong, 
heat deposition effects near the surface, where there is 
a sharp density gradient, become important. The extra 
heat can change the opacity of the matter by increasing 
the ionization fraction, resulting in stronger amplifica- 
tion of the pulsation and in turn accelerating the pre- 
dicted timescale from the first pulsation until mass ejec- 
tion. Concurrently, however, shock dissipation can con- 
vert the kinetic energy of the envelope into heat, while 
mass loss can also efficiently remove kinetic energy from 
the star. Heat deposition, shock dissipation, and mass 
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loss interact in complicated ways to form the full dy- 
namical picture of the envelope in the later phases. 

Similar ejections of the outer layers in models of lumi- 
nous Cepheid models, however, are known to be related 
to numerics rather than physics; see Smolec (2016). Re- 
gardless, we do not follow the code until this phase be- 
cause the timestep becomes prohibitively small (~ 100 
s). In particular, numerical difficulties arise in the 
Newton-Raphson iterations, during which the code fails 
to resolve the formation of convection zones around the 
shock front. Due to shock compression, the extra heat- 
ing also invalidates the equilibrium assumptions of the 
mixing length theory (Vitense 1953). To prevent non- 
physical convective behavior from developing near the 
surface, we set a ceiling for the relative convective ve- 
locity with respect to local velocity and the speed of 
sound in the simulations (see also the Appendix for re- 
lated setting in MESA). To limit the steepness of shock- 
waves and distribute them over multiple zones, explicit 
pulsation codes like RSP include either artificial viscos- 
ity or eddy viscosity terms (or both), but this is only 
effective up to certain L/M ratios (Stellingwerf 1975; 
Smolec 2016). 

Also at this stage, non-linear effects become dominant, 
causing sub-annual features to appear gradually on top 
of the linear pulsation. As observations of Betelgeuse 
do not show periodicities on sub-annual timescales, we 
do not consider this phase of pulsation further, though 
a study of non-linear pulsation with the dynamical cou- 
pling of opacity and ionization will be interesting future 
work. 

By comparing the general features of our hydrody- 
namical model with the pulsation patterns of Betelgeuse, 
it becomes clear that the quasi-annual variation is in- 
deed caused solely by the contraction and expansion of 
the star. It is interesting to note that in this linear os- 
cillation phase, we do not see any evidence of longer 
timescale variations, such as the 6-year and 35-year pe- 
riodicities. In fact, the hydrodynamic simulations never 
reproduce any of the observed variations besides the cur- 
rently presented quasi-annual pulsation, even when the 
initial mass is varied. This implies that these periodici- 
ties are driven by some mechanism outside the scope of 
what 1D hydrodynamic simulations can reproduce, i.e., 
not the &-mechanism. It would be interesting to conduct 
further dynamical studies on how the star relaxes when 
the opacity effects becomes important; however, work 
in this domain will require an algorithm to suppress the 
development of the &-mechanism so that the pulsation 
can be sustained without triggering excessive mass loss. 

There are similar works in the literature that focus on 
the pulsational features of massive stellar envelopes us- 


ing the stellar evolution code described in Langer et al. 
(1988). In particular, Heger et al. (1997) present the 
dynamical evolution of massive stars from 10-20 Mo 
and analyze their linear stability. In Figure 18, we plot 
the phase diagram of our canonical model's log49(L/ Lo) 
against logo Tef during the hydrodynamic evolution as 
a means of comparing directly with Figure 5 in Heger 
et al. (1997). In their work, the 11 Mo Red Supergiant 
model is followed for about 75 periods of oscillation, 
whereas ours capture the first 45 periods. Beyond the 
45*^. our models show numerical instability where the 
expansion and compression interact with the convective 
mixing zone, which largely suppresses the timestep and 
creates non-linear behaviour. 

Heger et al. (1997) show an approximately circular tra- 
jectory that spirals outwards from log,) Ter (K) & 3.52 
and logįo(L/Lo) = 4.90, whereas our model shows an 
elliptical trajectory, vacillating between high L and high 
Tæ on one side and low L and low Teg on the other. 
The outward spiraling in our work and theirs demon- 
strates that both stars are undergoing dynamical insta- 
bility with a growing amplitude. As expected, Heger 
et al. (1997)’s model has a lower period because it is a 
lower-mass model. This implies a more compact enve- 
lope, which allows all 75 periods of oscillations to happen 
within 30 years—this is only half the time of our model. 

Both models show a clockwise trajectory. Since the 
radius, temperature, and luminosity are related by the 
blackbody radiation formula L = 4mopR?TA, this 
means that when the stellar models resume their ini- 
tial luminosities, the models achieve a higher maximum 
Teg (i.e., smaller R) and a lower minimum Teg (i.e., 
larger R). These features suggest that Tog, L and 
R achieve their local extrema simultaneously in Heger 
et al. (1997)'s model, but in our case, this relationship 
is slightly lagged. As shown in Figure 13, our model 
approaches its local extrema with a non-zero velocity; 
thus, the stellar radius, which affects Tog, reaches its 
local maximum and minimum later than L. Our cal- 
culations therefore reproduce the phase lag between the 
luminosity and the velocity that has been observed in 
many other, smaller pulsators before (Castor 1968; Sz- 
abó et al. 2007). 

In Yoon & Cantiello (2010), the hydrodynami- 
cal features of a 20 Mo star with a luminosity of 
log;9(L/Lo) = 5.05 and temperature of Teg = 3198 
K are analyzed. Their stellar parameters are similar to 
ours, where our hydrodynamical model assumes a 21 
Mo star with a slightly lower initial hydrostatic lumi- 
nosity at logj9(D/ Lo) = 5.01 and Teg = 4000 K. They 
model about 50 years of the stellar pulsation; Figure 2 
in their work shows the surface velocity and is compa- 
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rable to Figure 13 in this work. Approximately 20 pul- 
sation cycles are followed in their work, where a higher 
period of ~ 1000 days is observed. Compared to our 
~ 400 day period, this indicates that their envelope is 
more relaxed and expanded. Both Yoon & Cantiello 
(2010) and our work show a consistent growth of the 
surface velocity. It takes about 5 cycles for the surface 
velocity to reach a ten-fold of amplification, while our 
model takes much longer—almost 20 cycles. This sug- 
gests that the &-mechanism is less efficient in our model, 
where the star exhibits behavior closer to adiabatic os- 
cillations than driven oscillations. From the growth of 
kinetic energy of the system, we can estimate that it 
takes a further ~ 10 years for the pulsation to grow 
to a surface velocity comparable with Yoon & Cantiello 
(2010). This would correspond to another 11-15 cy- 
cles in our case. Despite this, the robust exponential 
growth of the pulsation energy (see Figure 16) in both 
works implies that the pulsation could remove the out- 
ermost layers of the H-envelope from the star. However, 
this is not consistent with observational evidence; the 
RSGs we have observed pulsate with limited amplitude 
for several decades. The mass loss rate of the Betelgeuse 
is observed to be roughly ~ 1079 Mo yr! (Dolan et al. 
2016). We may speculate that driven pulsation is par- 
tially responsible for this driven mass loss, but the pul- 
sational growth rate is limited by the constant energy 
dissipation through shock heating and mass ejection. It 
is also true that the typical pulsation amplitude is not 
as large as would be expected for a driven wind (Fox & 
Wood 1982; Wood 2000; Wood et al. 2004). 

Whether or not mass loss can be driven depends on 
the degree of saturation in the pulsation of the surface 
layers (King et al. 1966). When the mode is permitted 
to develop, this process can be influential in the forma- 
tion of circumstellar matter in Type-IIn supernovae for 
massive stars close to ~ 20 Mo (e.g., Smith 2017). How- 
ever, given the regulated oscillation amplitude observed 
in a number of RSGs empirically, additional mechanisms 
not modeled in this work must become dominant in reg- 
ulating the growth of these oscillation patterns. 

'The most recent analysis of this kind can be found in 
Goldberg et al. (2020). The pulsation of a red super- 
giant with 16.3 Mo is computed using MESA with the 
GYRE extension (version 11701). In contrast to the ap- 
proaches discussed above, their hydrodynamic models 
involve an initial perturbation to the density distribu- 
tion to trigger direct pulsation of not only the funda- 
mental mode, but also the first overtone. They obtain 
a star of log;9(L/ Lo) = 5.2 and 880R5, which is about 
1096 larger than our model. As a result, their pulsation 
shows a fundamental mode with a longer period—about 


600 days—and first overtone of ~ 300 days. We reiterate 
that since we do not perturb the density profile at the 
onset of hydrodynamic evolution, and we assume that 
all pulsation is triggered by numerical perturbations, it 
is consistent that we do not see higher order excited 
modes alongside the fundamental mode. However, the 
strength of the quasi-annual variation of Betelgeuse but 
absence of clear, shorter-timescale periods suggests that 
we should not be concerned by a lack of overtone activity 
in our hydrodynamic modeling. The results in Goldberg 
et al. (2020) demonstrate that an overtone in our model 
would give rise to significant sub-annual motion, which 
is not observed in Betelgeuse. Furthermore, the growth 
rate depends in part on the L/M ratio, and this fac- 
tor may not be high enough in our models to excite an 
overtone. 

Goldberg et al. (2020)'s density perturbation approach 
can also cause the star to pulsate with a much larger 
amplitude initially, which is not achieved in our work 
nor in the two analyses presented above (Heger et al. 
1997; Yoon & Cantiello 2010). Therefore, it is unclear if 
their work shows a similar exponential growth as in the 
literature and in our model. 


5.5. Impact of Initial Mass om Pulsation 


''hus far, we have presented one model with an initial 
mass of 21 Mo. Now, we consider a series of hydrody- 
namic models of different initial mass and discuss how 
the progenitor mass affects the pulsation pattern. 

We repeat the simulations by varying the progenitor 
mass, while fixing the mixing length parameter (see Sec- 
tion A for details on the exact configuration) so that we 
can compare consistently among models. In Table 3, 
we tabulate the global parameters and pulsation statis- 
tics of these models. The data present the following 
trends: when the progenitor mass increases, the present 
day mass Msn, helium core mass Mye, the radius at the 
end of He-burning R, and its corresponding luminosity 
L all increase. There is a weak decreasing trend for 
the effective temperature Tẹ. Meanwhile, the time re- 
quired for the non-linear pulsation to emerge decreases. 
We note a severe drop between models of 20.2 and 20.5 
Mo during which the associated number of pulsations 
also decreases sharply. Below M = 19 Mo, the oscilla- 
tion does not amplify significantly within 300 years, at 
which point we terminate the simulation. 

We note in particular two entries in Table 3 showing 
simulations with initial (evolutionary) masses of 20M: 
one with ayyr=3.0 and one with aypr=2.5. In the 
latter case, the non-linear excitation time drops consid- 
erably. This data point disrupts an otherwise monoton- 
ically decreasing trend in pulse number with increasing 
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Table 3. T'he global properties and pulsation statistics of the 
hydrodynamical models studied in this work. M, Ma, and 
Muye are the initial, final and He-core masses of the star in 
units of Mo, respectively. R, log, ) L and Teg are the initial 
radius in units of Ro, luminosity in units of Lo and the 
effective temperature in units of K at the beginning of the 
hydrodynamical phase, respectively. a is the mixing length 
parameter. trun is the time the pulsation of the star becomes 
non-linear, where we stop the simulations, in years. “Pulse” is 
the number of pulsation cycles experienced by the star before 
the onset of non-linear pulsations. No number is available 
when trun is larger than 300 years. All hydrodynamic models 
are launched from the evolutionary point at which the helium 
mass fraction in the core is 10 ^. Though this does not 
correspond to the evolutionary phase statistically preferred 
by our classical models, this starting condition lends itself to 
more stable hydrodynamic models and allows us to explore 
the appropriate radius, luminosity, and mass regimes. 


M Ma, Mye R logiyLl Teg lian Pulse 
18 17.12 5.57 550 4.92 4160 >300 N/A 
19 17.90 6.06 624 5.00 4115 >300 N/A 
20 18.80 6.47 655 5.04 4117 166.2 ~ 230 

18.95 6.60 659 5.05 4120 141.8 ~ 200 


19.17 6.75 707 5.10 4081 31.5 43 
19.54 7.00 721 5.11 4083 27.0 40 
20.30 7.46 787 5.18 4053 21.0 24 
20.92 8.04 875 5.25 4008 16.7 18 


w 
o 
ho 
wwww ww w w]e 


19 2.5 17.78 60.07 724 5.01 3832 102.8 122 


20 2.5 18.57 6.59 794 5.08 3801 41.1 37 


present-day mass (above 19.Mc), demonstrating that the 
relationship between pulse number, trun, and present- 
day mass is not totally straightforward. It is well-known, 
however, that below a certain level of precision, the im- 
pact on global parameters caused by changes in the 
mixing length are indistinguishable from variations in 
mass and metallicity in models of low to intermediate 
mass stars with convective envelopes (Joyce & Chaboyer 
2018a). Further, the late stage evolution of high mass 
stars is especially sensitive to convective parameters; in 
practice, ayr is often tuned arbitrarily until the model 
converges or behaves as desired. We include the last 
row of Table 3 to highlight this degeneracy and caution 
against over-interpretation. 

In Figure 17, we present the time evolution of the 
pulsation pattern for models with the progenitor mass 
from 19 to 22 Mo. We choose these masses as their 
timescales are more relevant to that of Betelgeuse. Be- 
fore non-linearity disturbs the simulation, all models be- 
have similarly in both luminosity and radius. 

Despite the fact that the excitation time apparently 
depends on the progenitor mass and mixing length pa- 


rameter, the the means by which the star becomes 
excited—i.e. the pulsation driving mechanism itself—is 
less sensitive to these choices. The dynamical pulsation 
always concludes with a significant drop in the stellar 
luminosity, and the peak luminosity and maximum ra- 
dius are similar among all models near the end of the 
simulation. 

'To further outline the similarity, we plot in Figure 18 
the phase diagram of representative models from 19- 
23 Mo. Clear similarity can be seen for models above 
19 Mo. In particular, for M = 20 Mo, the model 
has a highly extended trun of ~ 166 year. All models 
have an elliptical structure, which is actually a clock- 
wise outward-going spiral. They show once again that 
all stars evolve toward a high L and a high Teg state si- 
multaneously, or the converse. This suggests that the 
driving mechanism in all of these models is qualita- 
tively the same, too. A higher progenitor mass gives 
rise to a sparser trajectory; however, we notice that for 
M = 19 Mo, there is no regularity in the trajectory. 
This suggests that the &-mechanism fails to stimulate 
residual numerical noise into periodic motions. 

'To further characterize the runaway time of the M — 
20 Mo model, we compare the amount of time needed 
for the star to develop non-linear pulsation (tun) after 
using the hydrodynamical prescription. For progenitor 
masses above 20.5 Mo, trun decreases slowly with time. 
As shown in Figure 17, nonlinear activation timescales 
for masses of this range are between 15-40 years. How- 
ever, below 20.5 Mo, there is a sudden jump in trun, 
and the star requires more than roughly > 150 in or- 
der for non-linearity to become significant. The sudden 
jump could signify some qualitative changes in the stel- 
lar profile, namely that the «mechanism becomes much 
less effective in amplifying the acoustic wave inside the 
star; a detailed comparison to and analysis of the means 
of formation for the &-mechanism will be an interesting 
future project, but is beyond the scope of the present 
work. Crucially, this mass-sensitive timescale bifurca- 
tion suggests that the time required for the star to de- 
velop non-linear pulsation could be a highly discerning 
attribute among models of Betelgeuse. 

As an order of magnitude estimation, the typical lu- 
minosity of our model star is 10° Lo. The amount of 
energy dissipated is then ~ 1046 erg per year, but the 
kinetic energy is only on the order of 10*!74? erg. This is 
because radiation acts as a damping force through pho- 
ton emission, and without a consistent driving force for 
the pulsation, the oscillation would quickly dissipate. 

In the previous text, we have shown that there are 
multiple periodicities in Betelgeuse's lightcurve. These 
include a quasi-annual mode, a 6-year mode, a 30-year 
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Figure 17. The temporal evolution of the luminosity and stellar radius, scaled by their initial values, is shown for models with 
different progenitor masses studied in this work. Time 0 marks the transition from hydrostatic stellar evolutionary calculations 


to the hydrodynamical prescription. See also Appendix A for the exact numerical treatment. 


modulation, and, potentially, an overtone mode with 
185 day period. In the hydrodynamic models, we re- 
cover only the 416 day period. These results are largely 
self-consistent, as the 416 d mode is driven by the « 
mechanism, the LSP is not, and the 30 year modula- 
tion is most likely caused by rotation, which is not an 
internally driven form of variability. In the case of the 
overtone, however, we must address the question of how 
multiple modes may appear in the first place. 

One possibility is by non-linear mode excitation, as 
touched upon in Section 5.3. Through large amplitude 
oscillations, the outer layers can accumulate sufficient 
energy and momentum to compress matter beneath the 
stellar surface. This results in compression heating, 
which in turn raises the local temperature. This may im- 
pact the convective structure in the near-surface regions, 
thus presenting an additional source of energy that al- 
ters the net energy flow inside the star. Capturing this 
scenario numerically is particularly challenging because 
it involves modeling the dynamics of mixing behaviour 


in the convection zone. Meanwhile, the standard mixing 
length theory adopted in our work assumes the convec- 
tive mixing is in equilibrium (Vitense 1953). Modeling 
this phase properly would require a more sophisticated 
approach to time-dependent mixing and a robust solv- 
ing mechanism. We reiterate that the development of 
the overtone mode can also be sensitive to the numeri- 
cal setting. In particular, the implicit nature of the hy- 
drodynamics tends to suppress acoustic waves naturally, 
regardless of the use of artificial viscosity. "Therefore, 
the current non-detection can be attributed to numerics 
alone. Future exploration using, for example, an explicit 
hydrodynamical scheme would shed light on this matter 
but is beyond the scope of this study. 

We observe that in our hydrodynamical grid, models 
with M > 21 Mg develop non-linear pulsation much 
more quickly than the lower mass models, as shown in 
Figure 17. We emphasize that our results only suggest 
that the star is developing non-linear pulsation and/or 
near-surface shocks; how strong the final shock is re- 


UNCOVERING a ORIONIS 27 


T T 
M-19M,, 4 5.05 


sui 


log L/L. 
log, (L/L. 


| 5.04} E 
l 1 ji 1 j 1 l 1 l 1 
3.616 3.614 3.612 3.618 3.616 3.614 3.612 
logio (Tael K) log, (T al K) 
5.124 
51- | F l 
"i | e if i 
2 | a | | 
E | Ep 4 
E E | 
5.09} [ | 
l L l il 1 l il ] 
3.614 3.612 3.61 3.608 3.614 3.612 3.61 3.608 
log, [um IK) log, (Tay! K) 
; : 
5.18} M=22 Mon H [ MUS Mus | 
~ ~ 5.25} 4 
E 4| | 
a a I | 
$ $n | 
E E L ] 
l 1 l 1 ji L ji L ji L l 
ud 3.61 3.608 3.606 Jed 3.604 3.602 3.6 
log, Tae” K) logo (Lay! K) 


Figure 18. The phase diagrams for the model with a progenitor mass 19 (top left), 20 (top right), 20.5 (middle left), 21 (middle 
right), 22 (bottom left), 23 (bottom right) Mo respectively. The trajectory is cut when the non-linearity begins to disturb the 
elliptical pattern in each figure. 
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mains unclear. Because we do not follow our hydrody- 
namical simulations far enough to investigate rigorously 
the various energy dissipation channels discussed in sec- 
tion 5.4, it remains to be studied whether o Orionis’s 
strong pulsations eventually reach a new equilibrium or 
behave in some other way. 

Another possible excitation mechanism is wave-driven 
pulsation, as described in Shiode & Quataert (2014); 
Fuller (2017); Fuller & Ro (2018). This mechanism 
proposes that a convective wave in the star can par- 
tially penetrate through the evanescent regions! and 
approach the stellar surface. Although wave-driven pul- 
sation was described in the context of very late phases 
of stellar evolution in those works (i.e.; Neon-Oxygen 
burning, rather than He), the theory suggests that as 
long as convection is activated, energy can be trans- 
ferred from the interior convection zone to regions near 
the surface, where it can then excite surface motion. 
However, depending on the convective luminosity, such 
a mechanism would provide a heavily condensed energy 
deposition near the surface, in turn triggering enormous 
losses in mass of 0.01-1 Mo yr-!. Mass loss of this order 
is not observed in Betelgeuse. 


6. CONCLUSIONS 


We have presented a detailed observational and theo- 
retical analysis of o Orionis, including the presentation 
of new photometry and three different types of numer- 
ical predictions from classical evolutionary, linear seis- 
mic, and hydrodynamic simulations. Our critical results 
are summarized as follows. 

We present a new set of processed, space-based pho- 
tometric data from the SMEI instrument, filling a gap 
in precise, publicly available photometry during the 
late 2000s. These data reveal variation on monthly 
timescales, which is likely the signature of convective 
cell turnover. In combination with longitudinal data 
collected by the AAVSO, the photometry confirms the 
presence of several key periodicities and contextualizes 
the recent dimming behavior of o Orionis in the long- 
term. 

We determine that the fundamental mode and the 
LSP are longer according to the SMEI and ground-based 
V-band photometric data than according to the long- 
term visual results of Kiss et al. (2006): Po = 416 +24 d 
and Psp = 2365+10 d in our work versus P = 388-:30 
d and Ppsp ~ 2050 d in theirs. We conclude that the 
semiregular variability of the star—except the primary 
dimming event of 2019-2020—can be explained by phase 
changes in a short-lifetime pulsation mode and the pho- 
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tometric effects of giant convective cells. We also detect 
a new component with 185 + 13.5 d period. We iden- 
tify this as the first overtone, thus classifying o Ori a 
double-mode pulsator. 

We conduct a grid-based analysis of evolutionary 
tracks to estimate the fundamental, model-derived pa- 
rameters of œ Ori. Supported by previous studies, 
we take special account of the theoretical uncertainty 
imparted by an ad hoc choice of the mixing length 
parameter, awpr, and reconsider the uncertainties on 
Betelgeuse’s effective temperature accordingly (Joyce & 
Chaboyer 2018a,b; Levesque & Massey 2020). We per- 
form a probabilistic age prior analysis and find good 
agreement between our estimates of Betelgeuse’s cur- 
rent evolutionary stage (RSB core helium burning) and 
present-day mass range (16.5-19 Mo) with previous 
modeling initiatives (Neilson et al. 2011; Dolan et al. 
2016; Wheeler et al. 2017; Nance et al. 2018). Our 
seismic analysis prefers a median initial mass range of 
18—21Mọo. However, we find that the observed, present- 
day rotational velocity of a Ori cannot be reproduced us- 
ing single-star evolution; a merger or some other source 
of spin-up is required, in agreement with Wheeler et al. 
(2017); Chatzopoulos et al. (2020). The likelihood of a 
previous interaction is also supported by our kinematic 
argument in Section 2. 

Linear seismic analysis with GYRE heavily constrains 
the radius of Betelgeuse, for which we report a value of 
764-8 Ro. Combining this result with existing angu- 
lar diameter and temperature data, we obtain a parallax 
value of 7 = 5.95*025 mas for a Orionis based on seis- 


mic constraints, resulting in a precise and independent 
distance estimate of 168*72 pc. Our results are con- 
sistent with reprocessed Hipparcos measurements but 
in disagreement with recent radio parallax observations 
(van Leeuwen 2007; Harper et al. 2017), highlighting 
the difficulty of estimating cosmic noise when deriving 
the geometric parallax of this star. To the best of our 
knowledge, this is the first time that a seismic parallax 
has been obtained for Betelgeuse. 

Deeper analysis of emergent periodicities in both hy- 
drostatic seismic and hydrodynamic models, in conjunc- 
tion with existing observational data on variable stars 
across the mass spectrum, unambiguously demonstrate 
that the 416 d period derived in this work is due to 
pulsation in the fundamental p-mode. 

Finally, using hydrodynamic models with six differ- 
ent masses, we investigate the physics of these oscilla- 
tions. All hydrodynamic models in the prescribed mass 
range manifest similar quasi-annual behavior as the fun- 
damental mode, in agreement with similar studies. Our 
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hydrodynamical simulations thus confirm that the 416 
day pulsation is driven by the &-mechanism. 

We find that stars with an initial mass below ~ 20 Mo 
take much longer for the pulsation to excite other oscil- 
lation modes; in particular, a 19 Mo model can take 
as long as 150 years to build up to non-linearity. The 
similarity among models suggests that the exact param- 
eters of the model play a less important role in repro- 
ducing the fundamental mode of the star. Importantly, 
if non-linear excitation is assumed to be correlated to 
the &-mechanism's triggering of overtone modes, and if 
the observed mass loss in Betelgeuse is not pulsationally 
driven, our hydrodynamic simulations constrain against 
progenitor masses above ~ 20 Mo. On the contrary, 
if the large amplitude pulsation fails to reach an equi- 
librium and instead triggers shock waves and consecu- 
tive mass loss, it would strongly suggest that Betelgeuse 
has a mass greater than 19M. As we do not evolve 
our simulations far enough to characterize the late-stage 
pulsational behavior, we cannot infer mass constraints 
definitively from these simulations. 

It is unclear whether the excited fundamental mode 
can be modulated by other radiative mechanisms or lead 
to observable mass loss. If mass loss can be triggered, 
the short runaway time from the appearance of the first 
wave until mass ejection suggests that the star can lose 
a considerable amount of its H-envelope during its post- 
main-sequence evolution. If the observed mass loss in 
Betelgeuse can be connected to the instability observed 
in this work, we could potentially make additional infer- 
ences about the initial mass of Betelgeuse based on the 
timescale of non-linear excitation. 

The sudden bifurcation in excitation time as a 
function of mass in our hydrodynamical models pro- 
vides some constraint on Betelgeuse's upcoming, pre- 
supernova evolution. For models with an initial mass 
above ~ 20 Mg (present-day mass 18.8 Mo), the &- 
mechanism driven pulsation and the mass loss it incites 
could partially remove the H-envelope prior to the final 
explosion. This would give rise to a Type-IIp, Type-IIL 
and then Type-IIn supernova. Meanwhile, for models 
with initial masses below this break-off point, the very 
long excitation time of the k-mechanism means that the 
star would retain most of its H-envelope. In this case, 
an alternative mass loss channel would be required for 
the formation of a circumstellar medium. 

Conclusively determining which of these two possi- 
ble evolutionary channels o Ori will take would require 
disentangling the degeneracy between mass and mixing 
length in the simulations, but our work here suggests 
that a predictive investigation in this vein is possible. 
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APPENDIX 


A. MESA CONFIGURATIONS 


In this section, we detail the configuration profile for the evolutionary and hydrodynamical portions of the simulations. 
The evolutionary phase inherits settings from the massive_star_defaults inlist. Additionally, we set the “Dutch” 
mass loss prescription with a parameter 0.8, namely: 


hot_wind_scheme = ’Dutch’ 
Dutch_scaling factor = 0.8 
hot wind full on T = 1d0 
hot wind full off T = OdO 


In order to construct a star that maintains the proper radius for hydrodynamic evolution, we must adjust the mixing 
length parameter: 


mixing length alpha - 3 
MLT. opion = ’Henyey’ 


We notice that a larger mixing length parameter results in a smaller radius at the end of the He-burning. The mass of 
the star is selected such that the luminosity is within the expected range (~ 4.8-5.1) and a radius between 700-800 R5 
for consistency with the seismic parameters. 

A requirement of our configuration is that the star should exhibit an observable amount of pulsation within a 
reasonable amount of time (~ 100 years). A small progenitor mass results in a very long quiescent time. Meanwhile, 
a higher mass can trigger observable pulsation quickly, but its luminosity and radius can be too high. As a result, 
for the hydrodynamics, we pick the high mass end M = 21 Mọ with a large mixing length parameter a = 3. This is 
slightly higher than what is used in the evolutionary calculations (a < 2.5), but our model gives the correct radius at 
720 Ro and a luminosity ~ 109-4 Lo. The final mass (present-day mass) is 19.5 Mo and the helium core is 7.00Mo 

As we require that the stellar profile transition smoothly from the evolutionary phase to the hydrodynamical phase, 
we use identical settings in the dynamical phase. 


T mix limit - O 
min. T for acceleration limited, conv velocity = 0 
okay. to reduce. gradT excess = .false. 


In the hydrodynamical phase, we patch extra settings onto this configuration such that the hydrostatic equilibrium 
constructed in the previous phase can be well maintained. However, one qualitative change is included, where the 
mass loss is suspended. 


Dutch scaling factor = 0.0d0 


This is a reasonable approximation given that we are simulating a short period of time: ~ 100 years. 
To trigger the hydrodynamics, we use the standard settings as provided by the test, suite test case ccsn in MESA 
version 8118. This includes 


use. ODE var.egn pairing = .true. 

use dvdt form of momentum eqn = .true. 

use. dPrad, dm, form of T gradient eqn = .true. 
use dedt form of energy .eqgn = .true. 

use momentum, outer BC - .true. 

use ODE form of density .eqgn = .true. 


These settings have been used in our previous work modeling the dynamical pulsation in pulsation pair-instability 
supernovae. See Leung et al. (2019, 2020) for the application of these setting to the more massive star counterpart. 

Furthermore, to ensure the code captures the early oscillation when the simulation has begun, we impose a maximum 
evolutionary timestep of 10? s. 
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max_timestep = 100000 


We also remove the temperature limitation in which the hydrodynamics is solved. This means the Euler equations are 
solved throughout the star, without assuming the envelope is in hydrostatic equilibrium: 


velocity_logT_lower_bound = 0 


To prevent supersonic convection from occurring in the simulation and invalidating the assumptions of the mixing 


length theory, we impose a cap on the convective speed via 


mlt_accel_g_ theta = 1 
max_v_div_cs_for_convection = 1.0d-1 
max_conv_vel_div_csound = 1.0d0, 


ensuring that the convective behavior remains physical. 


At last, we turn on the artificial viscosity so that all potential shocks can be resolved by the simulation. This 
happens, in particular, near the surface where the density gradient is the highest. 


use_artificial_viscosity = .true. 
shock_spread_linear = 0 
shock_spread_quadratic = 2d-2 


We find that a higher artificial viscosity parameter can result in the code crashing earlier in the simulation, whereas a 
value too small can result in too strong of a shock when the global pulsation amplitude is still weak. 
A simulation of ~ 30 years requires approximately 10000 timesteps. 
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